Fixed parallel partitioning and a create leaf bug

common/sys/vector_t.h

@@ -23,22 +23,22 @@ namespace embree
 
 #if defined(VECTOR_INIT_ALLOCATOR)
     template<typename M>
-    vector_t (M alloc) 
+    vector_t (M alloc)
       : alloc(alloc), size_active(0), size_alloced(0), items(nullptr) {}
 
     template<typename M>
-    vector_t (M alloc, size_t sz) 
+    vector_t (M alloc, size_t sz)
       : alloc(alloc), size_active(0), size_alloced(0), items(nullptr) { resize(sz); }
 
 #else
-      vector_t () 
+      vector_t ()
         : size_active(0), size_alloced(0), items(nullptr) {}
-    
-      vector_t (size_t sz) 
+
+      vector_t (size_t sz)
         : size_active(0), size_alloced(0), items(nullptr) { resize(sz); }
 #endif
 
-      
+
       ~vector_t() {
         clear();
       }
@@ -51,7 +51,7 @@ namespace embree
         for (size_t i=0; i<size_active; i++) items[i] = other.items[i];
       }
 
-      vector_t& operator=(const vector_t& other) 
+      vector_t& operator=(const vector_t& other)
       {
         resize(other.size_active);
         for (size_t i=0; i<size_active; i++) items[i] = other.items[i];
@@ -75,10 +75,10 @@ namespace embree
         internal_resize(new_size,size_alloced < new_size ? new_size : size_alloced);
       }
 
-      void reserve(size_t new_alloced) 
+      void reserve(size_t new_alloced)
       {
         /* do nothing if container already large enough */
-        if (new_alloced <= size_alloced) 
+        if (new_alloced <= size_alloced)
           return;
 
         /* resize exact otherwise */
@@ -103,10 +103,10 @@ namespace embree
       __forceinline       T* data()       { return items; };
       __forceinline const T* data() const { return items; };
 
-     
+
       /******************** Modifiers **************************/
 
-      __forceinline void push_back(const T& nt) 
+      __forceinline void push_back(const T& nt)
       {
         const T v = nt; // need local copy as input reference could point to this vector
         internal_grow(size_active+1);
@@ -118,12 +118,12 @@ namespace embree
         size_active--;
       }
 
-      void clear() 
+      void clear()
       {
         /* destroy elements */
         for (size_t i=0; i<size_active; i++)
           alloc.destroy(&items[i]);
-        
+
         /* free memory */
         alloc.deallocate(items,size_alloced); items = nullptr;
         size_active = size_alloced = 0;
@@ -133,7 +133,7 @@ namespace embree
 
       void internal_resize(size_t new_active, size_t new_alloced)
       {
-        assert(new_active <= new_alloced); 
+        assert(new_active <= new_alloced);
 
         /* destroy elements */
         for (size_t i=new_active; i<size_active; i++)
@@ -143,9 +143,9 @@ namespace embree
         size_active = new_active;
 
         /* only reallocate if necessary */
-        if (new_alloced == size_alloced) 
+        if (new_alloced == size_alloced)
           return;
-        
+
         /* reallocate and copy items */
         T* old_items = items;
         items = alloc.allocate(new_alloced);
@@ -157,7 +157,7 @@ namespace embree
       void internal_grow(size_t new_alloced)
       {
         /* do nothing if container already large enough */
-        if (new_alloced <= size_alloced) 
+        if (new_alloced <= size_alloced)
           return;
 
         /* resize to next power of 2 otherwise */

kernels/algorithms/IndexSort.cpp

@@ -0,0 +1,366 @@
+//
+//  IndexSort.cpp
+//  RayAccelerator
+//
+//  Created by Rasmus Barringer on 2014-09-03.
+//  Copyright (c) 2014 Rasmus Barringer. All rights reserved.
+//
+
+#include "IndexSort.h"
+#include <immintrin.h>
+#include <algorithm>
+
+void sortedIndicesFromFloats(float* __restrict values, unsigned* __restrict indices, unsigned count) {
+	long long* keys = reinterpret_cast<long long*>(indices);
+
+	__m256i i0 = _mm256_setr_epi32(0, 1, 2, 3, 4, 5, 6, 7);
+	__m256i i1 = _mm256_setr_epi32(8, 9, 10, 11, 12, 13, 14, 15);
+	__m256i i2 = _mm256_setr_epi32(16, 17, 18, 19, 20, 21, 22, 23);
+	__m256i i3 = _mm256_setr_epi32(24, 25, 26, 27, 28, 29, 30, 31);
+
+	__m256i increment = _mm256_set1_epi32(32);
+	__m256i signMask = _mm256_set1_epi32(0x80000000);
+
+	unsigned i = 0;
+
+	for (; i+31 < count; i += 32) {
+		const __m256i* input = reinterpret_cast<const __m256i*>(values + i);
+		__m256i* output = reinterpret_cast<__m256i*>(keys + i);
+
+		__m256i value0 = _mm256_load_si256(input + 0);
+		__m256i value1 = _mm256_load_si256(input + 1);
+		__m256i value2 = _mm256_load_si256(input + 2);
+		__m256i value3 = _mm256_load_si256(input + 3);
+
+		__m256i absValue0 = _mm256_andnot_si256(signMask, value0);
+		__m256i absValue1 = _mm256_andnot_si256(signMask, value1);
+		__m256i absValue2 = _mm256_andnot_si256(signMask, value2);
+		__m256i absValue3 = _mm256_andnot_si256(signMask, value3);
+
+		__m256i encodedValue0 = _mm256_sign_epi32(absValue0, value0);
+		__m256i encodedValue1 = _mm256_sign_epi32(absValue1, value1);
+		__m256i encodedValue2 = _mm256_sign_epi32(absValue2, value2);
+		__m256i encodedValue3 = _mm256_sign_epi32(absValue3, value3);
+
+		__m256i key00 = _mm256_unpacklo_epi32(i0, encodedValue0);
+		__m256i key01 = _mm256_unpackhi_epi32(i0, encodedValue0);
+		__m256i key10 = _mm256_unpacklo_epi32(i1, encodedValue1);
+		__m256i key11 = _mm256_unpackhi_epi32(i1, encodedValue1);
+		__m256i key20 = _mm256_unpacklo_epi32(i2, encodedValue2);
+		__m256i key21 = _mm256_unpackhi_epi32(i2, encodedValue2);
+		__m256i key30 = _mm256_unpacklo_epi32(i3, encodedValue3);
+		__m256i key31 = _mm256_unpackhi_epi32(i3, encodedValue3);
+
+		_mm256_store_si256(output + 0, key00);
+		_mm256_store_si256(output + 1, key01);
+		_mm256_store_si256(output + 2, key10);
+		_mm256_store_si256(output + 3, key11);
+		_mm256_store_si256(output + 4, key20);
+		_mm256_store_si256(output + 5, key21);
+		_mm256_store_si256(output + 6, key30);
+		_mm256_store_si256(output + 7, key31);
+
+		i0 = _mm256_add_epi32(i0, increment);
+		i1 = _mm256_add_epi32(i1, increment);
+		i2 = _mm256_add_epi32(i2, increment);
+		i3 = _mm256_add_epi32(i3, increment);
+	}
+
+	for (; i < count; ++i) {
+		int input = *reinterpret_cast<const int*>(values + i);
+		int absInput = input & 0x7fffffff;
+
+		if (input < 0)
+			absInput = -absInput;
+
+		keys[i] = ((long long)absInput << 32) | i;
+	}
+
+	std::sort(keys, keys + count);
+
+	for (unsigned i = 0; i < count; i += 32) {
+		const float* input = reinterpret_cast<const float*>(keys + i);
+		float* output = reinterpret_cast<float*>(indices + i);
+
+		__m256 value0 = _mm256_load_ps(input + 0);
+		__m256 value1 = _mm256_load_ps(input + 8);
+		__m256 value2 = _mm256_load_ps(input + 16);
+		__m256 value3 = _mm256_load_ps(input + 24);
+		__m256 value4 = _mm256_load_ps(input + 32);
+		__m256 value5 = _mm256_load_ps(input + 40);
+		__m256 value6 = _mm256_load_ps(input + 48);
+		__m256 value7 = _mm256_load_ps(input + 56);
+
+		__m256 i0 = _mm256_shuffle_ps(value0, value1, _MM_SHUFFLE(2,0,2,0));
+		__m256 i1 = _mm256_shuffle_ps(value2, value3, _MM_SHUFFLE(2,0,2,0));
+		__m256 i2 = _mm256_shuffle_ps(value4, value5, _MM_SHUFFLE(2,0,2,0));
+		__m256 i3 = _mm256_shuffle_ps(value6, value7, _MM_SHUFFLE(2,0,2,0));
+
+		i0 = _mm256_castpd_ps(_mm256_permute4x64_pd(_mm256_castps_pd(i0), (0) | ((2) << 2) | ((1) << 4) | ((3) << 6)));
+		i1 = _mm256_castpd_ps(_mm256_permute4x64_pd(_mm256_castps_pd(i1), (0) | ((2) << 2) | ((1) << 4) | ((3) << 6)));
+		i2 = _mm256_castpd_ps(_mm256_permute4x64_pd(_mm256_castps_pd(i2), (0) | ((2) << 2) | ((1) << 4) | ((3) << 6)));
+		i3 = _mm256_castpd_ps(_mm256_permute4x64_pd(_mm256_castps_pd(i3), (0) | ((2) << 2) | ((1) << 4) | ((3) << 6)));
+
+		_mm256_store_ps(output + 0, i0);
+		_mm256_store_ps(output + 8, i1);
+		_mm256_store_ps(output + 16, i2);
+		_mm256_store_ps(output + 24, i3);
+	}
+}
+
+inline void radixSortUint64High32X3(unsigned long long* __restrict dataA, unsigned long long* __restrict dataB, unsigned long long* __restrict dataC,
+									unsigned long long* __restrict tempA, unsigned long long* __restrict tempB, unsigned long long* __restrict tempC, unsigned count) {
+	unsigned bucketsA0[257];
+	unsigned bucketsA1[257];
+	unsigned bucketsA2[257];
+	unsigned bucketsA3[257];
+
+	unsigned bucketsB0[257];
+	unsigned bucketsB1[257];
+	unsigned bucketsB2[257];
+	unsigned bucketsB3[257];
+
+	unsigned bucketsC0[257];
+	unsigned bucketsC1[257];
+	unsigned bucketsC2[257];
+	unsigned bucketsC3[257];
+
+	for (unsigned i = 0; i < 257; ++i) {
+		bucketsA0[i] = 0;
+		bucketsA1[i] = 0;
+		bucketsA2[i] = 0;
+		bucketsA3[i] = 0;
+
+		bucketsB0[i] = 0;
+		bucketsB1[i] = 0;
+		bucketsB2[i] = 0;
+		bucketsB3[i] = 0;
+
+		bucketsC0[i] = 0;
+		bucketsC1[i] = 0;
+		bucketsC2[i] = 0;
+		bucketsC3[i] = 0;
+	}
+
+	unsigned* histogramA0 = bucketsA0+1;
+	unsigned* histogramA1 = bucketsA1+1;
+	unsigned* histogramA2 = bucketsA2+1;
+	unsigned* histogramA3 = bucketsA3+1;
+
+	unsigned* histogramB0 = bucketsB0+1;
+	unsigned* histogramB1 = bucketsB1+1;
+	unsigned* histogramB2 = bucketsB2+1;
+	unsigned* histogramB3 = bucketsB3+1;
+
+	unsigned* histogramC0 = bucketsC0+1;
+	unsigned* histogramC1 = bucketsC1+1;
+	unsigned* histogramC2 = bucketsC2+1;
+	unsigned* histogramC3 = bucketsC3+1;
+
+	for (unsigned i = 0; i < count; ++i) {
+		unsigned long long da = dataA[i];
+		unsigned long long db = dataB[i];
+		unsigned long long dc = dataC[i];
+
+		++histogramA0[(da >> (4 << 3)) & 0xff];
+		++histogramA1[(da >> (5 << 3)) & 0xff];
+		++histogramA2[(da >> (6 << 3)) & 0xff];
+		++histogramA3[(da >> (7 << 3)) & 0xff];
+
+		++histogramB0[(db >> (4 << 3)) & 0xff];
+		++histogramB1[(db >> (5 << 3)) & 0xff];
+		++histogramB2[(db >> (6 << 3)) & 0xff];
+		++histogramB3[(db >> (7 << 3)) & 0xff];
+
+		++histogramC0[(dc >> (4 << 3)) & 0xff];
+		++histogramC1[(dc >> (5 << 3)) & 0xff];
+		++histogramC2[(dc >> (6 << 3)) & 0xff];
+		++histogramC3[(dc >> (7 << 3)) & 0xff];
+	}
+
+	for (unsigned i = 1; i < 256; ++i) {
+		bucketsA0[i] += bucketsA0[i-1];
+		bucketsA1[i] += bucketsA1[i-1];
+		bucketsA2[i] += bucketsA2[i-1];
+		bucketsA3[i] += bucketsA3[i-1];
+
+		bucketsB0[i] += bucketsB0[i-1];
+		bucketsB1[i] += bucketsB1[i-1];
+		bucketsB2[i] += bucketsB2[i-1];
+		bucketsB3[i] += bucketsB3[i-1];
+
+		bucketsC0[i] += bucketsC0[i-1];
+		bucketsC1[i] += bucketsC1[i-1];
+		bucketsC2[i] += bucketsC2[i-1];
+		bucketsC3[i] += bucketsC3[i-1];
+	}
+
+	for (unsigned i = 0; i < count; ++i) {
+		unsigned long long da = dataA[i];
+		unsigned long long db = dataB[i];
+		unsigned long long dc = dataC[i];
+
+		unsigned indexA = bucketsA0[(da >> (4 << 3)) & 0xff]++;
+		unsigned indexB = bucketsB0[(db >> (4 << 3)) & 0xff]++;
+		unsigned indexC = bucketsC0[(dc >> (4 << 3)) & 0xff]++;
+
+		tempA[indexA] = da;
+		tempB[indexB] = db;
+		tempC[indexC] = dc;
+	}
+
+	for (unsigned i = 0; i < count; ++i) {
+		unsigned long long da = tempA[i];
+		unsigned long long db = tempB[i];
+		unsigned long long dc = tempC[i];
+
+		unsigned indexA = bucketsA1[(da >> (5 << 3)) & 0xff]++;
+		unsigned indexB = bucketsB1[(db >> (5 << 3)) & 0xff]++;
+		unsigned indexC = bucketsC1[(dc >> (5 << 3)) & 0xff]++;
+
+		dataA[indexA] = da;
+		dataB[indexB] = db;
+		dataC[indexC] = dc;
+	}
+
+	for (unsigned i = 0; i < count; ++i) {
+		unsigned long long da = dataA[i];
+		unsigned long long db = dataB[i];
+		unsigned long long dc = dataC[i];
+
+		unsigned indexA = bucketsA2[(da >> (6 << 3)) & 0xff]++;
+		unsigned indexB = bucketsB2[(db >> (6 << 3)) & 0xff]++;
+		unsigned indexC = bucketsC2[(dc >> (6 << 3)) & 0xff]++;
+
+		tempA[indexA] = da;
+		tempB[indexB] = db;
+		tempC[indexC] = dc;
+	}
+
+	for (unsigned i = 0; i < count; ++i) {
+		unsigned long long da = tempA[i];
+		unsigned long long db = tempB[i];
+		unsigned long long dc = tempC[i];
+
+		unsigned indexA = bucketsA3[(da >> (7 << 3)) & 0xff]++;
+		unsigned indexB = bucketsB3[(db >> (7 << 3)) & 0xff]++;
+		unsigned indexC = bucketsC3[(dc >> (7 << 3)) & 0xff]++;
+
+		dataA[indexA] = da;
+		dataB[indexB] = db;
+		dataC[indexC] = dc;
+	}
+}
+
+inline void packFloatsAndIndicesInUint64(float* __restrict values, unsigned long long* __restrict keys, unsigned count, int offset) {
+
+	__m256i offset256 = _mm256_set1_epi32(offset);
+	__m256i i0 = _mm256_add_epi32(_mm256_setr_epi32(0, 1, 2, 3, 4, 5, 6, 7), offset256);
+	__m256i i1 = _mm256_add_epi32(_mm256_setr_epi32(8, 9, 10, 11, 12, 13, 14, 15), offset256);
+	__m256i i2 = _mm256_add_epi32(_mm256_setr_epi32(16, 17, 18, 19, 20, 21, 22, 23), offset256);
+	__m256i i3 = _mm256_add_epi32(_mm256_setr_epi32(24, 25, 26, 27, 28, 29, 30, 31), offset256);
+
+	__m256i increment = _mm256_set1_epi32(32);
+	__m256 signMask = _mm256_castps_si256(_mm256_set1_epi32(0x80000000));
+	__m256 fullMask = _mm256_castps_si256(_mm256_set1_epi32(0xffffffff));
+
+	unsigned i = 0;
+
+	for (; i+31 < count; i += 32) {
+		const float* input = reinterpret_cast<const float*>(values + i);
+		float* output = reinterpret_cast<float*>(keys + i);
+
+		__m256 value0 = _mm256_load_ps(input + 0*8);
+		__m256 value1 = _mm256_load_ps(input + 1*8);
+		__m256 value2 = _mm256_load_ps(input + 2*8);
+		__m256 value3 = _mm256_load_ps(input + 3*8);
+
+		__m256 mask0 = _mm256_blendv_ps(signMask, fullMask, value0);
+		__m256 mask1 = _mm256_blendv_ps(signMask, fullMask, value1);
+		__m256 mask2 = _mm256_blendv_ps(signMask, fullMask, value2);
+		__m256 mask3 = _mm256_blendv_ps(signMask, fullMask, value3);
+
+		__m256 encodedValue0 = _mm256_xor_ps(mask0, value0);
+		__m256 encodedValue1 = _mm256_xor_ps(mask1, value1);
+		__m256 encodedValue2 = _mm256_xor_ps(mask2, value2);
+		__m256 encodedValue3 = _mm256_xor_ps(mask3, value3);
+
+		__m256 key00 = _mm256_unpacklo_ps(i0, encodedValue0);
+		__m256 key01 = _mm256_unpackhi_ps(i0, encodedValue0);
+		__m256 key10 = _mm256_unpacklo_ps(i1, encodedValue1);
+		__m256 key11 = _mm256_unpackhi_ps(i1, encodedValue1);
+		__m256 key20 = _mm256_unpacklo_ps(i2, encodedValue2);
+		__m256 key21 = _mm256_unpackhi_ps(i2, encodedValue2);
+		__m256 key30 = _mm256_unpacklo_ps(i3, encodedValue3);
+		__m256 key31 = _mm256_unpackhi_ps(i3, encodedValue3);
+
+		_mm256_store_ps(output + 0*8, key00);
+		_mm256_store_ps(output + 1*8, key01);
+		_mm256_store_ps(output + 2*8, key10);
+		_mm256_store_ps(output + 3*8, key11);
+		_mm256_store_ps(output + 4*8, key20);
+		_mm256_store_ps(output + 5*8, key21);
+		_mm256_store_ps(output + 6*8, key30);
+		_mm256_store_ps(output + 7*8, key31);
+
+		i0 = _mm256_add_epi32(i0, increment);
+		i1 = _mm256_add_epi32(i1, increment);
+		i2 = _mm256_add_epi32(i2, increment);
+		i3 = _mm256_add_epi32(i3, increment);
+	}
+
+	for (; i < count; ++i) {
+		int input = *reinterpret_cast<const int*>(values + i);
+
+		input ^= input < 0 ? 0xffffffff : 0x80000000;
+
+		keys[i] = ((unsigned long long)input << 32) | (i + offset);
+	}
+}
+
+inline void inPlaceUnpackIndicesFromLowUin64ToUint32(unsigned long long* keys, unsigned count) {
+	for (unsigned i = 0; i < count; i += 32) {
+		const float* input = reinterpret_cast<const float*>(keys + i);
+		float* output = reinterpret_cast<float*>(keys) + i;
+
+		__m256 value0 = _mm256_load_ps(input + 0);
+		__m256 value1 = _mm256_load_ps(input + 8);
+		__m256 value2 = _mm256_load_ps(input + 16);
+		__m256 value3 = _mm256_load_ps(input + 24);
+		__m256 value4 = _mm256_load_ps(input + 32);
+		__m256 value5 = _mm256_load_ps(input + 40);
+		__m256 value6 = _mm256_load_ps(input + 48);
+		__m256 value7 = _mm256_load_ps(input + 56);
+
+		__m256 i0 = _mm256_shuffle_ps(value0, value1, _MM_SHUFFLE(2,0,2,0));
+		__m256 i1 = _mm256_shuffle_ps(value2, value3, _MM_SHUFFLE(2,0,2,0));
+		__m256 i2 = _mm256_shuffle_ps(value4, value5, _MM_SHUFFLE(2,0,2,0));
+		__m256 i3 = _mm256_shuffle_ps(value6, value7, _MM_SHUFFLE(2,0,2,0));
+
+		i0 = _mm256_permute4x64_pd(i0, (0) | ((2) << 2) | ((1) << 4) | ((3) << 6));
+		i1 = _mm256_permute4x64_pd(i1, (0) | ((2) << 2) | ((1) << 4) | ((3) << 6));
+		i2 = _mm256_permute4x64_pd(i2, (0) | ((2) << 2) | ((1) << 4) | ((3) << 6));
+		i3 = _mm256_permute4x64_pd(i3, (0) | ((2) << 2) | ((1) << 4) | ((3) << 6));
+
+		_mm256_store_ps(output + 0, i0);
+		_mm256_store_ps(output + 8, i1);
+		_mm256_store_ps(output + 16, i2);
+		_mm256_store_ps(output + 24, i3);
+	}
+}
+
+void sortedIndicesFromFloatsX3(float* __restrict valuesA, float* __restrict valuesB, float* __restrict valuesC, unsigned* __restrict indicesA, unsigned* __restrict indicesB, unsigned* __restrict indicesC, unsigned count, int offset) {
+	unsigned long long* keysA = reinterpret_cast<unsigned long long*>(indicesA);
+	unsigned long long* keysB = reinterpret_cast<unsigned long long*>(indicesB);
+	unsigned long long* keysC = reinterpret_cast<unsigned long long*>(indicesC);
+
+	packFloatsAndIndicesInUint64(valuesA, keysA, count, offset);
+	packFloatsAndIndicesInUint64(valuesB, keysB, count, offset);
+	packFloatsAndIndicesInUint64(valuesC, keysC, count, offset);
+
+	radixSortUint64High32X3(keysA, keysB, keysC, keysA + count, keysB + count, keysC + count, count);
+
+	inPlaceUnpackIndicesFromLowUin64ToUint32(keysA, count);
+	inPlaceUnpackIndicesFromLowUin64ToUint32(keysB, count);
+	inPlaceUnpackIndicesFromLowUin64ToUint32(keysC, count);
+}

kernels/algorithms/IndexSort.h

@@ -0,0 +1,24 @@
+//
+//  IndexSort.h
+//  RayAccelerator
+//
+//  Created by Rasmus Barringer on 2014-09-03.
+//  Copyright (c) 2014 Rasmus Barringer. All rights reserved.
+//
+
+#ifndef RayAccelerator_IndexSort_h
+#define RayAccelerator_IndexSort_h
+
+// Note: indices needs to be 2x in size for temporary storage.
+// Everything needs to be aligned and padded:
+//  values - pad with 128 bytes
+//  indices - pad with 256 bytes
+void sortedIndicesFromFloats(float* __restrict values, unsigned* __restrict indices, unsigned count);
+
+// Note: indices needs to be 4x in size for temporary storage.
+// Everything needs to be aligned and padded:
+//  values - pad with 128 bytes
+//  indices - pad with 256 bytes
+void sortedIndicesFromFloatsX3(float* __restrict valuesA, float* __restrict valuesB, float* __restrict valuesC, unsigned* __restrict indicesA, unsigned* __restrict indicesB, unsigned* __restrict indicesC, unsigned count, int offset);
+
+#endif

kernels/common/scene.cpp

@@ -24,7 +24,7 @@
 #include "../xeonphi/bvh4mb/bvh4mb.h"
 #include "../xeonphi/bvh4hair/bvh4hair.h"
 #endif
- 
+
 namespace embree
 {
   /* error raising rtcIntersect and rtcOccluded functions */
@@ -34,26 +34,26 @@ namespace embree
   void invalid_rtcIntersect8() { throw_RTCError(RTC_INVALID_OPERATION,"rtcIntersect8 and rtcOccluded8 not enabled"); }
   void invalid_rtcIntersect16() { throw_RTCError(RTC_INVALID_OPERATION,"rtcIntersect16 and rtcOccluded16 not enabled"); }
 
- 
+
 /* number of created scene */
   AtomicCounter Scene::numScenes = 0;
 
   Scene::Scene (Device* device, RTCSceneFlags sflags, RTCAlgorithmFlags aflags)
-    : device(device), 
+    : device(device),
       Accel(AccelData::TY_UNKNOWN),
-      flags(sflags), aflags(aflags), numMappedBuffers(0), is_build(false), modified(true), 
-      needTriangleIndices(false), needTriangleVertices(false), 
+      flags(sflags), aflags(aflags), numMappedBuffers(0), is_build(false), modified(true),
+      needTriangleIndices(false), needTriangleVertices(false),
       needBezierIndices(false), needBezierVertices(false),
       needSubdivIndices(false), needSubdivVertices(false),
-      numTriangles(0), numTriangles2(0), 
-      numBezierCurves(0), numBezierCurves2(0), 
-      numSubdivPatches(0), numSubdivPatches2(0), 
+      numTriangles(0), numTriangles2(0),
+      numBezierCurves(0), numBezierCurves2(0),
+      numSubdivPatches(0), numSubdivPatches2(0),
       numUserGeometries1(0), numSubdivEnableDisableEvents(0),
       numIntersectionFilters4(0), numIntersectionFilters8(0), numIntersectionFilters16(0),
-      commitCounter(0), commitCounterSubdiv(0), 
+      commitCounter(0), commitCounterSubdiv(0),
       progress_monitor_function(nullptr), progress_monitor_ptr(nullptr), progress_monitor_counter(0)
   {
-#if defined(TASKING_LOCKSTEP) 
+#if defined(TASKING_LOCKSTEP)
     lockstep_scheduler.taskBarrier.init(MAX_MIC_THREADS);
 #elif defined(TASKING_TBB_INTERNAL)
     scheduler = nullptr;
@@ -88,16 +88,16 @@ namespace embree
     {
       std::cout << "scene flags: static " << isStatic() << " compact = " << isCompact() << " high quality = " << isHighQuality() << " robust = " << isRobust() << std::endl;
     }
-    
-    if (device->tri_accel == "default" || device->tri_accel == "bvh4i")   
+
+    if (device->tri_accel == "default" || device->tri_accel == "bvh4i")
     {
-      if (device->tri_builder == "default") 
+      if (device->tri_builder == "default")
       {
         if (isStatic())
         {
           if (device->verbosity(1)) std::cout << "STATIC BUILDER MODE" << std::endl;
           if ( isCompact() )
-            accels.add(BVH4i::BVH4iTriangle1MemoryConservativeBinnedSAH(this,isRobust()));		    
+            accels.add(BVH4i::BVH4iTriangle1MemoryConservativeBinnedSAH(this,isRobust()));
           else if ( isHighQuality() )
             accels.add(BVH4i::BVH4iTriangle1ObjectSplitBinnedSAH(this,isRobust()));
           else
@@ -130,12 +130,12 @@ namespace embree
         else if (device->tri_builder == "morton64") {
           accels.add(BVH4i::BVH4iTriangle1ObjectSplitMorton64Bit(this,isRobust()));
         }
-        
+
         else THROW_RUNTIME_ERROR("unknown builder "+device->tri_builder+" for BVH4i<Triangle1>");
       }
     }
     else THROW_RUNTIME_ERROR("unknown accel "+device->tri_accel);
-    
+
 #else
     createTriangleAccel();
     accels.add(BVH4::BVH4Triangle4vMB(this));
@@ -153,32 +153,32 @@ namespace embree
 
   void Scene::createTriangleAccel()
   {
-    if (device->tri_accel == "default") 
+    if (device->tri_accel == "default")
     {
       if (isStatic()) {
-        int mode =  2*(int)isCompact() + 1*(int)isRobust(); 
+        int mode =  2*(int)isCompact() + 1*(int)isRobust();
         switch (mode) {
-        case /*0b00*/ 0: 
+        case /*0b00*/ 0:
 #if defined (__TARGET_AVX__)
           if (hasISA(AVX))
 	  {
             // FIXME: there appears to be a bug in the spatial split builder, performance and SAH is lower than for standard builder
-            //if (isHighQuality()) accels.add(BVH8::BVH8Triangle4SpatialSplit(this)); 
-            /*else*/               accels.add(BVH8::BVH8Triangle4ObjectSplit(this)); 
+            //if (isHighQuality()) accels.add(BVH8::BVH8Triangle4SpatialSplit(this));
+            /*else*/               accels.add(BVH8::BVH8Triangle4ObjectSplit(this));
           }
-          else 
+          else
 #endif
           {
             // FIXME: there appears to be a bug in the spatial split builder, performance and SAH is lower than for standard builder
-            //if (isHighQuality()) accels.add(BVH4::BVH4Triangle4SpatialSplit(this));           
-            /*else*/               accels.add(BVH4::BVH4Triangle4ObjectSplit(this)); 
+            //if (isHighQuality()) accels.add(BVH4::BVH4Triangle4SpatialSplit(this));
+            /*else*/               accels.add(BVH4::BVH4Triangle4ObjectSplit(this));
           }
           break;
 
-          case /*0b01*/ 1: 
+          case /*0b01*/ 1:
 //#if defined (__TARGET_AVX2__) && !defined(__WIN32__) // FIXME: have to disable under Windows as watertightness tests fail
 //            if (hasISA(AVX2))
-//              accels.add(BVH8::BVH8Triangle8vObjectSplit(this)); 
+//              accels.add(BVH8::BVH8Triangle8vObjectSplit(this));
 //            else
 //#endif
               accels.add(BVH4::BVH4Triangle4vObjectSplit(this));
@@ -187,8 +187,8 @@ namespace embree
         case /*0b10*/ 2: accels.add(BVH4::BVH4Triangle4iObjectSplit(this)); break;
         case /*0b11*/ 3: accels.add(BVH4::BVH4Triangle4iObjectSplit(this)); break;
         }
-      } 
-      else 
+      }
+      else
       {
         int mode =  2*(int)isCompact() + 1*(int)isRobust();
         switch (mode) {
@@ -209,6 +209,7 @@ namespace embree
     else if (device->tri_accel == "bvh4.bvh4.triangle8")    accels.add(BVH4::BVH4BVH4Triangle8ObjectSplit(this));
     else if (device->tri_accel == "bvh4.triangle8")         accels.add(BVH4::BVH4Triangle8(this));
     else if (device->tri_accel == "bvh8.triangle4")         accels.add(BVH8::BVH8Triangle4(this));
+    else if (device->tri_accel == "bvh8.triangle4.bonsai")  accels.add(BVH8::BVH8Triangle4Bonsai(this));
     else if (device->tri_accel == "bvh8.triangle8")         accels.add(BVH8::BVH8Triangle8(this));
     else if (device->tri_accel == "bvh8.trianglepairs8")    accels.add(BVH8::BVH8TrianglePairs8(this));
     //else if (device->tri_accel == "bvh8.triangle8v")    accels.add(BVH8::BVH8Triangle8v(this));
@@ -219,18 +220,18 @@ namespace embree
 
   void Scene::createHairAccel()
   {
-    if (device->hair_accel == "default") 
+    if (device->hair_accel == "default")
     {
       if (isStatic()) {
-        int mode =  2*(int)isCompact() + 1*(int)isRobust(); 
+        int mode =  2*(int)isCompact() + 1*(int)isRobust();
         switch (mode) {
         case /*0b00*/ 0: accels.add(BVH4::BVH4OBBBezier1v(this,isHighQuality())); break;
         case /*0b01*/ 1: accels.add(BVH4::BVH4OBBBezier1v(this,isHighQuality())); break;
         case /*0b10*/ 2: accels.add(BVH4::BVH4OBBBezier1i(this,isHighQuality())); break;
         case /*0b11*/ 3: accels.add(BVH4::BVH4OBBBezier1i(this,isHighQuality())); break;
         }
-      } 
-      else 
+      }
+      else
       {
         int mode =  2*(int)isCompact() + 1*(int)isRobust();
         switch (mode) {
@@ -239,7 +240,7 @@ namespace embree
         case /*0b10*/ 2: accels.add(BVH4::BVH4Bezier1i(this)); break;
         case /*0b11*/ 3: accels.add(BVH4::BVH4Bezier1i(this)); break;
         }
-      }   
+      }
     }
     else if (device->hair_accel == "bvh4.bezier1v"    ) accels.add(BVH4::BVH4Bezier1v(this));
     else if (device->hair_accel == "bvh4.bezier1i"    ) accels.add(BVH4::BVH4Bezier1i(this));
@@ -250,7 +251,7 @@ namespace embree
 
   void Scene::createSubdivAccel()
   {
-    if (device->subdiv_accel == "default") 
+    if (device->subdiv_accel == "default")
     {
       if (isIncoherent(flags) && isStatic())
         accels.add(BVH4::BVH4SubdivGridEager(this));
@@ -265,7 +266,7 @@ namespace embree
 
 #endif
 
-  Scene::~Scene () 
+  Scene::~Scene ()
   {
     for (size_t i=0; i<geometries.size(); i++)
       delete geometries[i];
@@ -281,18 +282,18 @@ namespace embree
   void Scene::clear() {
   }
 
-  unsigned Scene::newUserGeometry (size_t items) 
+  unsigned Scene::newUserGeometry (size_t items)
   {
     Geometry* geom = new UserGeometry(this,items);
     return geom->id;
   }
-  
+
   unsigned Scene::newInstance (Scene* scene) {
     Geometry* geom = new Instance(this,scene);
     return geom->id;
   }
 
-  unsigned Scene::newTriangleMesh (RTCGeometryFlags gflags, size_t numTriangles, size_t numVertices, size_t numTimeSteps) 
+  unsigned Scene::newTriangleMesh (RTCGeometryFlags gflags, size_t numTriangles, size_t numVertices, size_t numTimeSteps)
   {
     if (isStatic() && (gflags != RTC_GEOMETRY_STATIC)) {
       throw_RTCError(RTC_INVALID_OPERATION,"static scenes can only contain static geometries");
@@ -303,12 +304,12 @@ namespace embree
       throw_RTCError(RTC_INVALID_OPERATION,"only 1 or 2 time steps supported");
       return -1;
     }
-    
+
     Geometry* geom = new TriangleMesh(this,gflags,numTriangles,numVertices,numTimeSteps);
     return geom->id;
   }
 
-  unsigned Scene::newSubdivisionMesh (RTCGeometryFlags gflags, size_t numFaces, size_t numEdges, size_t numVertices, size_t numEdgeCreases, size_t numVertexCreases, size_t numHoles, size_t numTimeSteps) 
+  unsigned Scene::newSubdivisionMesh (RTCGeometryFlags gflags, size_t numFaces, size_t numEdges, size_t numVertices, size_t numEdgeCreases, size_t numVertexCreases, size_t numHoles, size_t numTimeSteps)
   {
     if (isStatic() && (gflags != RTC_GEOMETRY_STATIC)) {
       throw_RTCError(RTC_INVALID_OPERATION,"static scenes can only contain static geometries");
@@ -324,14 +325,14 @@ namespace embree
 #if defined(__TARGET_AVX__)
     if (hasISA(AVX))
       geom = new SubdivMeshAVX(this,gflags,numFaces,numEdges,numVertices,numEdgeCreases,numVertexCreases,numHoles,numTimeSteps);
-    else 
+    else
 #endif
       geom = new SubdivMesh(this,gflags,numFaces,numEdges,numVertices,numEdgeCreases,numVertexCreases,numHoles,numTimeSteps);
     return geom->id;
   }
 
 
-  unsigned Scene::newBezierCurves (RTCGeometryFlags gflags, size_t numCurves, size_t numVertices, size_t numTimeSteps) 
+  unsigned Scene::newBezierCurves (RTCGeometryFlags gflags, size_t numCurves, size_t numVertices, size_t numTimeSteps)
   {
     if (isStatic() && (gflags != RTC_GEOMETRY_STATIC)) {
       throw_RTCError(RTC_INVALID_OPERATION,"static scenes can only contain static geometries");
@@ -342,17 +343,17 @@ namespace embree
       throw_RTCError(RTC_INVALID_OPERATION,"only 1 or 2 time steps supported");
       return -1;
     }
-    
+
     Geometry* geom = new BezierCurves(this,gflags,numCurves,numVertices,numTimeSteps);
     return geom->id;
   }
 
-  unsigned Scene::add(Geometry* geometry) 
+  unsigned Scene::add(Geometry* geometry)
   {
     Lock<AtomicMutex> lock(geometriesMutex);
 
     if (usedIDs.size()) {
-      int id = usedIDs.back(); 
+      int id = usedIDs.back();
       usedIDs.pop_back();
       geometries[id] = geometry;
       return id;
@@ -361,8 +362,8 @@ namespace embree
       return geometries.size()-1;
     }
   }
-  
-  void Scene::remove(Geometry* geometry) 
+
+  void Scene::remove(Geometry* geometry)
   {
     Lock<AtomicMutex> lock(geometriesMutex);
     usedIDs.push_back(geometry->id);
@@ -399,12 +400,12 @@ namespace embree
 
     /* select fast code path if no intersection filter is present */
     accels.select(numIntersectionFilters4,numIntersectionFilters8,numIntersectionFilters16);
-  
+
     /* build all hierarchies of this scene */
     accels.build(0,0);
-    
+
     /* make static geometry immutable */
-    if (isStatic()) 
+    if (isStatic())
     {
       accels.immutable();
       for (size_t i=0; i<geometries.size(); i++)
@@ -428,19 +429,19 @@ namespace embree
       std::cout << "selected scene intersector" << std::endl;
       intersectors.print(2);
     }
-    
+
     setModified(false);
   }
 
 #if defined(TASKING_LOCKSTEP)
 
-  void Scene::task_build_parallel(size_t threadIndex, size_t threadCount, size_t taskIndex, size_t taskCount, TaskScheduler::Event* event) 
+  void Scene::task_build_parallel(size_t threadIndex, size_t threadCount, size_t taskIndex, size_t taskCount, TaskScheduler::Event* event)
   {
     LockStepTaskScheduler::Init init(threadIndex,threadCount,&lockstep_scheduler);
     if (threadIndex == 0) accels.build(threadIndex,threadCount);
   }
 
-  void Scene::build (size_t threadIndex, size_t threadCount) 
+  void Scene::build (size_t threadIndex, size_t threadCount)
   {
     LockStepTaskScheduler::Init init(threadIndex,threadCount,&lockstep_scheduler);
     if (threadIndex != 0) return;
@@ -477,7 +478,7 @@ namespace embree
     }
 
     /* make static geometry immutable */
-    if (isStatic()) 
+    if (isStatic())
     {
       accels.immutable();
       for (size_t i=0; i<geometries.size(); i++)
@@ -508,13 +509,13 @@ namespace embree
 
 #if defined(TASKING_TBB_INTERNAL)
 
-  void Scene::build (size_t threadIndex, size_t threadCount) 
+  void Scene::build (size_t threadIndex, size_t threadCount)
   {
     AutoUnlock<MutexSys> buildLock(buildMutex);
 
     /* allocates own taskscheduler for each build */
     Ref<TaskSchedulerTBB> scheduler = nullptr;
-    { 
+    {
       Lock<MutexSys> lock(schedulerMutex);
       scheduler = this->scheduler;
       if (scheduler == null) {
@@ -553,7 +554,7 @@ namespace embree
 
 #if defined(TASKING_TBB)
 
-  void Scene::build (size_t threadIndex, size_t threadCount) 
+  void Scene::build (size_t threadIndex, size_t threadCount)
   {
     /* let threads wait for build to finish in rtcCommitThread mode */
     if (threadCount != 0) {
@@ -603,10 +604,10 @@ namespace embree
     unsigned int mxcsr = _mm_getcsr();
     _mm_setcsr(mxcsr | /* FTZ */ (1<<15) | /* DAZ */ (1<<6));
 #endif
-    
+
     try {
 
-#if TBB_INTERFACE_VERSION_MAJOR < 8    
+#if TBB_INTERFACE_VERSION_MAJOR < 8
       tbb::task_group_context ctx( tbb::task_group_context::isolated, tbb::task_group_context::default_traits);
 #else
       tbb::task_group_context ctx( tbb::task_group_context::isolated, tbb::task_group_context::default_traits | tbb::task_group_context::fp_settings );
@@ -622,14 +623,14 @@ namespace embree
           if (threadCount) group_barrier.wait(threadCount);
           group->wait();
 #if USE_TASK_ARENA
-        }); 
+        });
 #endif
-     
+
       /* reset MXCSR register again */
 #if !defined(__MIC__)
       _mm_setcsr(mxcsr);
 #endif
-    } 
+    }
     catch (...) {
 
        /* reset MXCSR register again */
@@ -656,7 +657,7 @@ namespace embree
     }
   }
 
-  void Scene::setProgressMonitorFunction(RTCProgressMonitorFunc func, void* ptr) 
+  void Scene::setProgressMonitorFunction(RTCProgressMonitorFunc func, void* ptr)
   {
     static MutexSys mutex;
     mutex.lock();

kernels/common/vector.h

@@ -35,17 +35,17 @@ namespace embree
       typedef std::size_t size_type; // FIXME: also use std::size_t type under windows if available
       typedef std::ptrdiff_t difference_type;
 
-      __forceinline aligned_monitored_allocator(MemoryMonitorInterface* device) 
+      __forceinline aligned_monitored_allocator(MemoryMonitorInterface* device)
         : device(device) {}
 
-      __forceinline pointer allocate( size_type n ) 
+      __forceinline pointer allocate( size_type n )
       {
         assert(device);
         device->memoryMonitor(n*sizeof(T),false);
         return (pointer) alignedMalloc(n*sizeof(value_type),alignment);
       }
 
-      __forceinline void deallocate( pointer p, size_type n ) 
+      __forceinline void deallocate( pointer p, size_type n )
       {
         assert(device);
         alignedFree(p);

kernels/xeon/CMakeLists.txt

@@ -53,12 +53,13 @@ IF (TASKING_TBB)
 ENDIF()
 
 SET(EMBREE_LIBRARY_FILES
- 
+
    embree.rc
 
   ../../common/tasking/taskscheduler_tbb.cpp
   ../../common/tasking/tasksys.cpp
 
+
   ../algorithms/parallel_for.cpp
   ../algorithms/parallel_reduce.cpp
   ../algorithms/parallel_prefix_sum.cpp
@@ -93,7 +94,7 @@ SET(EMBREE_LIBRARY_FILES
   geometry/instance_intersector1.cpp
   geometry/instance_intersector4.cpp
   geometry/subdivpatch1_intersector1.cpp
-  geometry/grid_soa.cpp		
+  geometry/grid_soa.cpp
   builders/primrefgen.cpp
 
   bvh4/bvh4.cpp
@@ -137,12 +138,12 @@ IF (RTCORE_RAY_PACKETS)
 ENDIF()
 
 SET(EMBREE_LIBRARY_FILES_AVX
-    
+
     ../common/scene_subdiv_mesh_avx.cpp
 
     geometry/primitive.cpp
     geometry/instance_intersector1.cpp
-    
+
     geometry/subdivpatch1_intersector1.cpp
     geometry/grid_soa.cpp
     builders/primrefgen.avx.cpp
@@ -155,10 +156,11 @@ SET(EMBREE_LIBRARY_FILES_AVX
     bvh4/bvh4_builder_twolevel.avx.cpp
     bvh4/bvh4_builder_subdiv.avx.cpp
     bvh4/bvh4_intersector1.cpp
-    
+
     bvh8/bvh8.cpp
     bvh8/bvh8_statistics.cpp
     bvh8/bvh8_builder_sah.avx.cpp
+
     bvh8/bvh8_intersector1.cpp
     )
 
@@ -182,13 +184,19 @@ ENDIF()
 
 SET(EMBREE_LIBRARY_FILES_AVX2
 
+
+
     geometry/instance_intersector1.cpp
     geometry/subdivpatch1_intersector1.cpp
     geometry/grid_soa.cpp
 
     bvh4/bvh4_intersector1.cpp
-    
-    bvh8/bvh8_intersector1.cpp)
+
+    bvh8/bvh8_intersector1.cpp
+builders/primrefgen.avx.cpp
+../algorithms/IndexSort.cpp
+bvh8/bvh8_builder_bonsai.cpp
+    )
 
 IF (RTCORE_RAY_PACKETS)
   SET(EMBREE_LIBRARY_FILES_AVX2 ${EMBREE_LIBRARY_FILES_AVX2}

kernels/xeon/builders/bvh_builder_bonsai.h

@@ -0,0 +1,908 @@
+// ======================================================================== //
+// Copyright 2009-2015 Intel Corporation                                    //
+//                                                                          //
+// Licensed under the Apache License, Version 2.0 (the "License");          //
+// you may not use this file except in compliance with the License.         //
+// You may obtain a copy of the License at                                  //
+//                                                                          //
+//     http://www.apache.org/licenses/LICENSE-2.0                           //
+//                                                                          //
+// Unless required by applicable law or agreed to in writing, software      //
+// distributed under the License is distributed on an "AS IS" BASIS,        //
+// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. //
+// See the License for the specific language governing permissions and      //
+// limitations under the License.                                           //
+// ======================================================================== //
+
+#pragma once
+
+#include "heuristic_bonsai.h"
+#include "../../common/scene.h"
+#include <chrono>
+
+namespace embree
+{
+  namespace isa
+  {
+    template<typename Set, typename Split>
+      struct GeneralBuildRecord
+      {
+      public:
+	__forceinline GeneralBuildRecord () {}
+
+        __forceinline GeneralBuildRecord (size_t depth)
+          : depth(depth), pinfo(empty) {}
+
+        __forceinline GeneralBuildRecord (const PrimInfo& pinfo, size_t depth, size_t* parent)
+          : pinfo(pinfo), depth(depth), parent(parent) {}
+
+        __forceinline GeneralBuildRecord (const PrimInfo& pinfo, size_t depth, size_t* parent, const Set &prims)
+          : pinfo(pinfo), depth(depth), parent(parent), prims(prims) {}
+
+          __forceinline GeneralBuildRecord (const PrimInfo& pinfo, size_t depth, size_t* parent, const Set &prims, const Set &indexRange)
+            : pinfo(pinfo), depth(depth), parent(parent), prims(prims), indexRange(indexRange) {}
+
+        __forceinline friend bool operator< (const GeneralBuildRecord& a, const GeneralBuildRecord& b) { return a.pinfo.size() < b.pinfo.size(); }
+	__forceinline friend bool operator> (const GeneralBuildRecord& a, const GeneralBuildRecord& b) { return a.pinfo.size() > b.pinfo.size(); }
+
+        __forceinline BBox3fa bounds() const { return pinfo.geomBounds; }
+
+        __forceinline size_t size() const { return this->pinfo.size(); }
+
+      public:
+        size_t* parent;   //!< Pointer to the parent node's reference to us
+		  size_t depth;     //!< Depth of the root of this subtree.
+		  Set prims;        //!< The list of primitives.
+		  PrimInfo pinfo;   //!< Bounding info of primitives.
+		  Split split;
+      Set indexRange;   //!< In memory available space for split partitioning
+      };
+
+    template<typename BuildRecord,
+      typename Heuristic,
+      typename ReductionTy,
+      typename Allocator,
+      typename CreateAllocFunc,
+      typename CreateNodeFunc,
+      typename UpdateNodeFunc,
+      typename CreateLeafFunc,
+      typename ProgressMonitor>
+
+      class GeneralBVHBuilder
+      {
+        static const size_t MAX_BRANCHING_FACTOR = 16;        //!< maximal supported BVH branching factor
+        static const size_t MIN_LARGE_LEAF_LEVELS = 8;        //!< create balanced tree of we are that many levels before the maximal tree depth
+        static const size_t SINGLE_THREADED_THRESHOLD = 4096; //!< threshold to switch to single threaded build
+		    size_t MAX_MINI_TREE_SIZE = 8192;		//!< maximum size of a mini tree
+		float  PRUNING_THRESHOLD = 0.95f;			//!< a pruning treshold in the range ]0,1[, anything else results in no pruning
+      const bool splitTriangles = false;
+      const int numberOfTriangles;
+		  volatile int miniTreeCount;
+      volatile int splitIndex;
+
+		  int underFilledNodes;
+      mvector<unsigned> ttIndicesX;
+      mvector<unsigned> ttIndicesY;
+      mvector<unsigned> ttIndicesZ;
+
+      mvector<float> ttMidpointsX;
+      mvector<float> ttMidpointsY;
+      mvector<float> ttMidpointsZ;
+
+      float* __restrict ttMidPoints[3];
+      unsigned* __restrict ttIndices[3];
+
+      mvector<PrimRef> ttBounds;
+
+      mvector<unsigned> ttTemporaryIndices;
+		  mvector<unsigned char> ttLeftPartition;
+		  mvector<float> ttAccumulatedArea;
+      mvector<unsigned> ttAccumulatedCounts;
+      mvector<FastSweepData*> fsd;
+      public:
+        mvector<BuildRecord> miniTreeRecords;
+		  mvector<BuildRecord> miniTreeRecordsPruned;
+      mvector<size_t*> miniTreeRoots;
+
+        GeneralBVHBuilder (Heuristic& heuristic,
+                           const ReductionTy& identity,
+                           CreateAllocFunc& createAlloc,
+                           CreateNodeFunc& createNode,
+                           UpdateNodeFunc& updateNode,
+                           CreateLeafFunc& createLeaf,
+                           ProgressMonitor& progressMonitor,
+                           const PrimInfo& pinfo,
+                           const size_t branchingFactor, const size_t maxDepth,
+                           const size_t logBlockSize, const size_t minLeafSize, const size_t maxLeafSize,
+                           const float travCost, const float intCost, const int numberOfTriangles, const bool splitTriangles, Scene* scene)
+          : heuristic(heuristic),
+          identity(identity),
+          createAlloc(createAlloc), createNode(createNode), updateNode(updateNode), createLeaf(createLeaf),
+          progressMonitor(progressMonitor),
+          pinfo(pinfo),
+          branchingFactor(branchingFactor), maxDepth(maxDepth),
+          logBlockSize(logBlockSize), minLeafSize(minLeafSize), maxLeafSize(maxLeafSize),
+          travCost(travCost), intCost(intCost), numberOfTriangles(numberOfTriangles),splitTriangles(splitTriangles),
+          ttTemporaryIndices(scene->device), ttLeftPartition(scene->device), ttAccumulatedArea(scene->device),
+          fsd(scene->device), miniTreeRecords(scene->device), miniTreeRecordsPruned(scene->device), miniTreeRoots(scene->device),
+          ttBounds(scene->device), ttIndicesX(scene->device), ttIndicesY(scene->device), ttIndicesZ(scene->device),
+          ttMidpointsX(scene->device), ttMidpointsY(scene->device), ttMidpointsZ(scene->device), ttAccumulatedCounts(scene->device)
+        {
+          if (branchingFactor > MAX_BRANCHING_FACTOR)
+            THROW_RUNTIME_ERROR("bvh_builder: branching factor too large");
+        }
+
+		  /*! bonsai builder */
+
+		  __forceinline void partition(BuildRecord& brecord, BuildRecord& lrecord, BuildRecord& rrecord) {
+			  heuristic.split(brecord.split,brecord.pinfo,brecord.prims,lrecord.pinfo,lrecord.prims,rrecord.pinfo, rrecord.prims);
+		  }
+
+
+
+		  __forceinline int getLargestDim(const BBox3fa& bounds) {
+			  Vec3fa dist = bounds.upper - bounds.lower;
+        int dim = -1;
+        if (dist[0] > 0 || dist[1] > 0 || dist[2] > 0) {
+			     dim = dist[0] > dist[1] ? 0 : 1;
+			      dim = dist[dim] > dist[2] ? dim : 2;
+          }
+
+			  return dim;
+		  }
+      __forceinline float getLargestDist(const BBox3fa& bounds) {
+			  Vec3fa dist = bounds.upper - bounds.lower;
+			    float largest = dist[0] > dist[1] ? dist[0] : dist[1];
+			      largest = largest > dist[2] ? largest : dist[2];
+
+
+			  return largest;
+		  }
+		  void miniTreeSelectRescursive(BuildRecord current) {
+
+
+        if (getLargestDist(current.pinfo.geomBounds) == 0)
+          return;
+			  if (current.prims.size() <= MAX_MINI_TREE_SIZE) {
+				  int miniTreeIndex = atomic_add(&miniTreeCount, 1);
+				  miniTreeRecords[miniTreeIndex] = current;
+				  return;
+			  }
+			  BuildRecord left, right;
+			  current.split.dim = getLargestDim(current.pinfo.centBounds);
+        current.split.pos = current.pinfo.centBounds.center()[current.split.dim];
+
+			  partition(current, left, right);
+
+
+			  if (current.size() > SINGLE_THREADED_THRESHOLD*4)
+			  {
+			  SPAWN_BEGIN;
+				SPAWN(([&] { miniTreeSelectRescursive(left); }));
+				SPAWN(([&] { miniTreeSelectRescursive(right); }));
+			  SPAWN_END;
+			  }
+			  else {
+				  miniTreeSelectRescursive(left);
+				  miniTreeSelectRescursive(right);
+			  }
+
+		  }
+
+		  __forceinline bool sweep(FastSweepData& data, BuildRecord& current, BuildRecord& left, BuildRecord& right, unsigned offset) {
+			  return heuristic.fastSweep(data, current.pinfo, current.split, current.prims, left.pinfo, left.prims, left.split, right.pinfo, right.prims, right.split, maxLeafSize, offset);
+		  }
+
+
+		  void buildMiniTreesRecursive(FastSweepData& data, BuildRecord& current, unsigned offset, Allocator alloc, size_t& miniTreeDepth) {
+
+			  if (alloc == nullptr)
+				 Allocator alloc = createAlloc();
+
+			  /*! initialize child list */
+			  BuildRecord children[MAX_BRANCHING_FACTOR];
+			  children[0] = current;
+			  size_t numChildren = 1;
+
+			  /*! split until node is full or SAH tells us to stop */
+			  do {
+
+				  /*! find best child to split */
+				  float largestSAH = -1;
+				  ssize_t largestChild = -1;
+				  for (size_t i=0; i<numChildren; i++)
+				  {
+
+					  if (children[i].pinfo.size() <= minLeafSize || children[i].split.isLeaf()) continue;
+					  float currentSAH = children[i].pinfo.leafSAH(logBlockSize);
+					  if (currentSAH > largestSAH) { largestChild = i; largestSAH = currentSAH; }
+
+				  }
+				  if (largestChild == -1) break;
+
+				  /* perform best found split */
+				  BuildRecord& brecord = children[largestChild];
+				  BuildRecord lrecord(current.depth+1);
+				  BuildRecord rrecord(current.depth+1);
+
+
+				  if (!sweep(data, brecord, lrecord, rrecord, offset)) {
+					  children[largestChild] = lrecord;
+					  children[numChildren] = rrecord;
+					  numChildren++;
+				  }
+
+			  } while (numChildren < branchingFactor);
+
+			  /*! create an inner node */
+        if (numChildren == 1) {
+          heuristic.gatherReferences(data, current.prims);
+          createLeaf(current, data.gatheredRefs, alloc);
+        return;
+        }
+        else
+          auto node = createNode(current,children,numChildren,alloc);
+
+        if (numChildren < 8) {
+          for (ssize_t i=numChildren-1; i>=0; i--) {
+
+				      heuristic.gatherReferences(data, children[i].prims);
+				      createLeaf(children[i], data.gatheredRefs, alloc);
+          }
+				  return;
+			  }
+			  /* spawn tasks */
+			  for (ssize_t i=numChildren-1; i>=0; i--) {
+          if (children[i].depth > miniTreeDepth)
+            miniTreeDepth = children[i].depth;
+				  buildMiniTreesRecursive(data, children[i], offset, alloc, miniTreeDepth);
+			  }
+		  }
+
+		  void setupIndices(unsigned** indices, const size_t first, const size_t count) {
+			  for (int i = 0; i < count; i++)
+				  indices[0][i] = indices[1][i] = indices[2][i] = i;// + first;
+		  }
+
+		  void buildMiniTrees(mvector<size_t*>& miniTreeRoots) {
+
+			  if (miniTreeCount == 1) {
+          FastSweepData data((int)(MAX_MINI_TREE_SIZE), maxLeafSize);
+          size_t depth = 0;
+  			  Allocator alloc = createAlloc();
+				  BuildRecord& br = miniTreeRecords[0];
+				  br.prims = range<size_t>(0, br.prims.size());
+				  heuristic.sortIndices(data, br.prims.size(), 0);
+				  buildMiniTreesRecursive(data, br, 0, alloc, depth);
+          br.depth = depth;
+			  }
+			  else {
+#if 0
+        FastSweepData data((int)(MAX_MINI_TREE_SIZE), maxLeafSize);
+        Allocator alloc = createAlloc();
+				  for (int i = 0; i < miniTreeCount; ++i) {
+					  BuildRecord& br = miniTreeRecords[i];
+					  data.setOffset(br.prims.begin());
+            size_t depth = 0;
+            unsigned offset = br.prims.begin();
+					  br.parent = (size_t*)&(miniTreeRoots.data()[i]);
+					  br.depth = 0;
+					  br.prims = range<size_t>(0, br.prims.size());
+					  br.pinfo.begin = 0;
+					  br.pinfo.end = br.prims.size();
+					  heuristic.sortIndices(data, br.prims.size(), offset);
+					  buildMiniTreesRecursive(data, br, offset, alloc, depth);
+            br.depth = depth;
+				  }
+#else
+
+          int threadCount = TaskSchedulerTBB::threadCount();
+
+          fsd.resize(threadCount);
+          FastSweepData** __restrict fastSweepDataPtr = fsd.data();
+          for (int i = 0; i < threadCount; i++) {
+              fastSweepDataPtr[i] = new FastSweepData(MAX_MINI_TREE_SIZE, maxLeafSize);
+          }
+
+  			  unsigned granularity = miniTreeCount/512;// / 128;
+
+          if (granularity < 1) {
+            granularity = miniTreeCount/128;
+          }
+          if (granularity < 1) {
+            granularity = 1;
+          }
+
+
+          size_t** __restrict miniTreeRootsPtr = miniTreeRoots.data();
+          //std::sort(miniTreeRecords.data(), miniTreeRecords.data() + miniTreeCount, [](BuildRecord& a, BuildRecord& b) { return a.prims.size() > b.prims.size(); });
+
+				  SPAWN_BEGIN;
+				  int treeID = 0;
+				  for (; treeID < miniTreeCount - granularity + 1; treeID += granularity) {
+
+            SPAWN(([this, treeID, miniTreeRootsPtr, fastSweepDataPtr, granularity] {
+              int threadIndex = TaskSchedulerTBB::threadIndex();
+              FastSweepData& data = *(fastSweepDataPtr[threadIndex]);
+
+              Allocator alloc = createAlloc();
+
+						  for (int i = treeID; i < treeID + granularity; ++i) {
+                BuildRecord& br = miniTreeRecords[i];
+                size_t depth = 0;
+                int offset = br.prims.begin();
+							  br.parent = (size_t*)&(miniTreeRootsPtr[i]);
+							  br.depth = 0;
+							  br.prims = range<size_t>(0, br.prims.size());
+							  br.pinfo.begin = 0;
+							  br.pinfo.end = br.prims.size();
+							  heuristic.sortIndices(data, br.prims.size(), offset);
+							  buildMiniTreesRecursive(data, br, offset, alloc, depth);
+                br.depth = depth;
+						}
+					  }));
+				  }
+          if (treeID < miniTreeCount) {
+				  SPAWN(([this, treeID, miniTreeRootsPtr, fastSweepDataPtr] {
+            int threadIndex = TaskSchedulerTBB::threadIndex();
+            FastSweepData& data = *(fastSweepDataPtr[threadIndex]);
+
+            Allocator alloc = createAlloc();
+
+					  for (int i = treeID; i < miniTreeCount; ++i) {
+						  BuildRecord& br = miniTreeRecords[i];
+						  //data.setOffset(br.prims.begin());
+              size_t depth = 0;
+              int offset = br.prims.begin();
+						  br.parent = (size_t*)&(miniTreeRootsPtr[i]);
+						  br.depth = 0;
+						  br.prims = range<size_t>(0, br.prims.size());
+						  br.pinfo.begin = 0;
+						  br.pinfo.end = br.prims.size();
+						  heuristic.sortIndices(data, br.prims.size(), offset);
+						  buildMiniTreesRecursive(data, br, offset, alloc, depth);
+              br.depth = depth;
+					  }}));
+          }
+        SPAWN_END;
+#endif
+			  }
+		}
+
+
+
+		  void partitionAccordingToRef(unsigned axis, unsigned first, unsigned last) {
+			  unsigned* __restrict temporaryIndices = ttTemporaryIndices.data();
+			  unsigned char* __restrict leftPartition = ttLeftPartition.data();
+
+			  //unsigned* __restrict indices = ttIndices[axis].data();
+        unsigned* __restrict indices = ttIndices[axis];
+
+			  unsigned* source = indices + first;
+			  unsigned* end = indices + last;
+
+			  unsigned* left = source;
+			  unsigned* right = temporaryIndices + first;
+
+			  while (source < end) {
+				  unsigned index = *source;
+				  unsigned l = leftPartition[index];
+
+				  *left = index;
+				  *right = index;
+
+				  left += l;
+				  right += 1-l;
+
+				  ++source;
+			  }
+
+			  std::copy(temporaryIndices + first, right, left);
+		  }
+
+		  void ttPartition(unsigned axis, unsigned first, unsigned last, unsigned pivot) {
+			  //unsigned* __restrict ref = (ttIndices[axis]).data();
+        unsigned* __restrict ref = (ttIndices[axis]);
+			  unsigned char* __restrict leftPartition = ttLeftPartition.data();
+
+			  for (unsigned i = first; i < pivot; ++i)
+				  leftPartition[ref[i]] = 1;
+
+			  for (unsigned i = pivot; i < last; ++i)
+				  leftPartition[ref[i]] = 0;
+
+			  partitionAccordingToRef((axis+1)%3, first, last);
+			  partitionAccordingToRef((axis+2)%3, first, last);
+		  }
+
+		  void topTreeSweep(BuildRecord& current, BuildRecord& lrecord, BuildRecord& rrecord) {
+        int block_shift = 3;
+			  unsigned first = current.prims.begin();
+			  unsigned last = current.prims.end();
+
+			  float bestSah = std::numeric_limits<float>::infinity();
+
+			  unsigned bestDim = 0xffffffff;
+			  unsigned pivot = 0xffffffff;
+
+			  BBox3fa bestBounds;
+
+        unsigned bestLeftCounts = 0;
+        unsigned bestRightCounts = 0;
+
+			  for (unsigned dim = 0; dim < 3; ++dim) {
+				  //unsigned* __restrict sortedIndices = ttIndices[dim].data();
+          unsigned* __restrict sortedIndices = ttIndices[dim];
+				  float* __restrict accumulatedArea = ttAccumulatedArea.data();
+          unsigned* __restrict accumulatedCounts = ttAccumulatedCounts.data();
+				  BBox3fa bounds = miniTreeRecordsPruned[sortedIndices[first]].pinfo.geomBounds;
+
+				  {
+					  int i = (int)first;
+            size_t lCounts = 0;
+					  for(; (int)i < (int)last-1; i++) {
+						  bounds.extend(miniTreeRecordsPruned[sortedIndices[i]].pinfo.geomBounds);
+						  float sa = halfArea(bounds);
+						  //lCounts += miniTreeRecordsPruned[sortedIndices[i]].pinfo.size();
+              lCounts = miniTreeRecordsPruned[sortedIndices[i]].depth > lCounts ? miniTreeRecordsPruned[sortedIndices[i]].depth : lCounts;
+              accumulatedArea[i] =  sa * sa * (float)((int)(i - first + 1));
+              //accumulatedArea[i] =  sa * float((size_t(i - first + 1)+(size_t(1)<<block_shift)-1) >> block_shift);
+              accumulatedCounts[i] = lCounts;
+					  }
+				  }
+
+				  bounds = miniTreeRecordsPruned[sortedIndices[last - 1]].pinfo.geomBounds;
+				  {
+					  int i = (int)last - 1;
+            size_t rCounts = 0;
+					  for (; (int)i > (int)first; i--) {
+						  bounds.extend(miniTreeRecordsPruned[sortedIndices[i]].pinfo.geomBounds);
+						  //rCounts += miniTreeRecordsPruned[sortedIndices[i]].pinfo.size();
+              rCounts = miniTreeRecordsPruned[sortedIndices[i]].depth > rCounts ? miniTreeRecordsPruned[sortedIndices[i]].depth : rCounts;
+              float sa = halfArea(bounds);
+						  float sah = (accumulatedArea[i - 1] + sa * sa * (float)((int)(last - i)));
+
+              //float sah = accumulatedArea[i - 1] +  sa * float((size_t(last-i)+(size_t(1)<<block_shift)-1) >> block_shift);
+
+						  if (sah < bestSah) {
+							  bestSah = sah;
+							  bestDim = dim;
+							  pivot = i;
+							  bestBounds = bounds;
+                bestLeftCounts = accumulatedCounts[i];
+                bestRightCounts = rCounts;
+						  }
+					  }
+				  }
+			  }
+
+			  lrecord.prims = range<size_t>(first, pivot);
+			  rrecord.prims = range<size_t>(pivot, last);
+
+			  ttPartition(bestDim, first, last, pivot);
+
+			  rrecord.pinfo = PrimInfo(pivot, last, bestBounds, bestBounds);
+        rrecord.depth = bestRightCounts;
+			  //bestBounds = miniTreeRecordsPruned[ttIndices[0].data()[first]].pinfo.geomBounds;
+        bestBounds = miniTreeRecordsPruned[ttIndices[0][first]].pinfo.geomBounds;
+			  size_t index = first + 1;
+			  for (; index < pivot; index++) {
+				  //bestBounds.extend(miniTreeRecordsPruned[ttIndices[0].data()[index]].pinfo.geomBounds);
+          bestBounds.extend(miniTreeRecordsPruned[ttIndices[0][index]].pinfo.geomBounds);
+			  }
+
+			  lrecord.pinfo = PrimInfo(first, pivot, bestBounds, bestBounds);
+        lrecord.depth = bestLeftCounts;
+		  }
+
+		  void buildTopTreeRecursive(BuildRecord& current, Allocator alloc) {
+
+			  if (alloc == nullptr)
+				  alloc = createAlloc();
+			  /* call memory monitor function to signal progress */
+			  //if (toplevel && current.size() <= SINGLE_THREADED_THRESHOLD)
+			  //	  progressMonitor(current.size());
+
+			  assert(current.pinfo.size() > 0);// || leafSAH >= 0 && splitSAH >= 0);
+
+        BuildRecord children[MAX_BRANCHING_FACTOR];
+			  /*! create a leaf node when threshold reached or SAH tells us to stop */
+			  if (current.pinfo.size() <= branchingFactor) {
+
+				  size_t count = current.pinfo.size();
+				  for (int i = 0; i < count; i++) {
+					   children[i] = miniTreeRecordsPruned[ttIndices[0][i + current.prims.begin()]];
+				  }
+
+				  if (count == 1) {
+					  *current.parent = *children[0].parent;
+
+				  }
+				  else {
+						  createNode(current, children, count, alloc, underFilledNodes);
+				  }
+
+				  return;
+			  }
+
+			  /*! initialize child list */
+
+			  children[0] = current;
+			  size_t numChildren = 1;
+
+			  /*! split until node is full or SAH tells us to stop */
+			  do {
+
+				  /*! find best child to split */
+				  float largestSAH = -1;
+				  ssize_t largestChild = -1;
+				  for (size_t i=0; i<numChildren; i++)
+				  {
+					  float currentSAH = children[i].pinfo.leafSAH(3);
+					  if (children[i].pinfo.size() == 1) continue;
+					  if (currentSAH > largestSAH) { largestChild = i; largestSAH = currentSAH; }
+				  }
+				  if (largestChild == -1) break;
+
+				  /* perform best found split */
+				  BuildRecord& brecord = children[largestChild];
+				  BuildRecord lrecord(current.depth+1);
+				  BuildRecord rrecord(current.depth+1);
+
+				  topTreeSweep(brecord, lrecord, rrecord);
+
+				  children[largestChild] = lrecord;
+				  children[numChildren] = rrecord;
+				  numChildren++;
+
+
+
+			  } while (numChildren < branchingFactor);
+
+			  /*! create an inner node */
+			  if (numChildren == 8)
+				  auto node = createNode(current,children,numChildren,alloc);
+
+			  /* spawn tasks */
+#if 0
+			  for (ssize_t i=numChildren-1; i>=0; i--) {
+				 // if (children[i].size() > 1)
+				  buildTopTreeRecursive(children[i], alloc);
+
+			  }
+#else
+			  if (current.size() > SINGLE_THREADED_THRESHOLD) {
+				  SPAWN_BEGIN;
+				  for (ssize_t i=numChildren-1; i>=0; i--)
+					  SPAWN(([&,i] { buildTopTreeRecursive(children[i],nullptr); }));
+				  SPAWN_END;
+			  }
+			  else {
+				  for (ssize_t i=numChildren-1; i>=0; i--) {
+					  buildTopTreeRecursive(children[i], alloc);
+
+				  }
+			  }
+#endif
+		  }
+
+      void sortTopTree() {
+
+        //if (!ttAlloc) {
+          Allocator alloc;
+          alloc = createAlloc();
+
+          ttIndicesX.resize(2*miniTreeCount+64);
+          ttIndicesY.resize(2*miniTreeCount+64);
+          ttIndicesZ.resize(2*miniTreeCount+64);
+
+          ttMidpointsX.resize(miniTreeCount+32);
+          ttMidpointsY.resize(miniTreeCount+32);
+          ttMidpointsZ.resize(miniTreeCount+32);
+
+          ttIndices[0] = (unsigned* __restrict)ttIndicesX.data();
+          ttIndices[1] = (unsigned* __restrict)ttIndicesY.data();
+          ttIndices[2] = (unsigned* __restrict)ttIndicesZ.data();
+
+        BuildRecord* __restrict mtrp = miniTreeRecordsPruned.data();
+        int i = 0;
+
+        float* __restrict mpx = ttMidpointsX.data();
+        float* __restrict mpy = ttMidpointsY.data();
+        float* __restrict mpz = ttMidpointsZ.data();
+
+        for (; i < miniTreeCount - 3; i += 4) {
+          BBox3fa b = mtrp[i].pinfo.geomBounds;
+				  __m128 mid0 = _mm_add_ps(b.lower, b.upper);
+          b = mtrp[i + 1].pinfo.geomBounds;
+					__m128 mid1 = _mm_add_ps(b.lower, b.upper);
+          b = mtrp[i + 2].pinfo.geomBounds;
+					__m128 mid2 = _mm_add_ps(b.lower, b.upper);
+          b = mtrp[i + 3].pinfo.geomBounds;
+					__m128 mid3 = _mm_add_ps(b.lower, b.upper);
+
+				  _MM_TRANSPOSE4_PS(mid0, mid1, mid2, mid3);
+				  _mm_store_ps(mpx + i, mid0);
+				  _mm_store_ps(mpy + i, mid1);
+				  _mm_store_ps(mpz + i, mid2);
+			  }
+
+			  for (; i < miniTreeCount; i++) {
+
+				  //Vec3fa mid = prims[(i + offset)].center2();
+          BBox3fa& b = mtrp[i].pinfo.geomBounds;
+					__m128 mid = _mm_add_ps(b.lower, b.upper);
+          ALIGNED(32) float midArray[4];
+					_mm_store_ps(midArray, mid);
+				  //prims[(i + offset)].lower = _mm_xor_ps(prims[(i + offset)].lower, neg);
+				  mpx[i] = midArray[0];
+				  mpy[i] = midArray[1];
+				  mpz[i] = midArray[2];
+
+			  }
+
+        SPAWN_BEGIN;
+        SPAWN(([&] {sortedIndicesFromFloats(mpx, ttIndices[0], miniTreeCount);}));
+        SPAWN(([&] {sortedIndicesFromFloats(mpy, ttIndices[1], miniTreeCount);}));
+        SPAWN(([&] {sortedIndicesFromFloats(mpz, ttIndices[2], miniTreeCount);}));
+        SPAWN_END;
+
+      }
+
+		  void buildTopTree(BuildRecord& rootRecord) {
+
+			  ttTemporaryIndices.resize(miniTreeCount);
+			  ttLeftPartition.resize(miniTreeCount);
+			  ttAccumulatedArea.resize(miniTreeCount);
+        ttAccumulatedCounts.resize(miniTreeCount);
+
+        sortTopTree();
+
+			  rootRecord.prims = range<size_t>(0, miniTreeCount);
+			  rootRecord.pinfo = PrimInfo(0, miniTreeCount, rootRecord.pinfo.geomBounds, rootRecord.pinfo.centBounds);
+			  rootRecord.depth = 0;
+        buildTopTreeRecursive(rootRecord, nullptr);
+		  }
+
+		  void pruneRecursive(int& newMiniTreeCount, const float averageArea, BVH8::NodeRef& nf, size_t depth) {
+			  depth -= 1;
+        for (int j = 0; j < branchingFactor; j++) {
+				  if (nf.node()->child(j).isNode() && (area(nf.node()->bounds(j)) > averageArea)) {
+					  pruneRecursive(newMiniTreeCount, averageArea, nf.node()->child(j), depth);
+				  }
+				  else if (nf.node()->child(j) != BVH8::emptyNode) {
+					  BuildRecord br;
+					  br.pinfo.geomBounds = nf.node()->bounds(j);
+					  br.parent = (size_t*)&nf.node()->child(j);
+            br.depth = depth;
+					  miniTreeRecordsPruned[newMiniTreeCount++] = br;
+				  }
+			  }
+		  }
+
+      // digs in to mini trees and removes too large pieces, which become mini trees of their own
+		  int prune() {
+			  float averageArea = 0;
+			  for (int i = 0; i < miniTreeCount; i++) {
+				  averageArea += area(miniTreeRecords[i].pinfo.geomBounds);
+
+			  }
+			  averageArea /= (float)miniTreeCount;
+
+        //std::cout << "area: " << averageArea << std::endl;
+			  averageArea *= PRUNING_THRESHOLD;
+			  int newMiniTreeCount = 0;
+			  int originalMiniTreeCount = miniTreeCount;
+        //std::cout << "originalMiniTreeCount: " << originalMiniTreeCount << std::endl;
+			  for (int i = 0; i < originalMiniTreeCount; i++) {
+          size_t originalDepth = miniTreeRecords[i].depth;
+
+				  BVH8::NodeRef nf = BVH8::NodeRef(*(miniTreeRecords[i].parent));
+				  if (nf.isNode() && area(miniTreeRecords[i].pinfo.geomBounds) > averageArea) {
+            originalDepth -= 1;
+					  for (int j = 0; j < branchingFactor; j++) {
+
+						  if (nf.node()->child(j).isNode() && (area(nf.node()->bounds(j)) > averageArea)) {
+							  pruneRecursive(newMiniTreeCount, averageArea, nf.node()->child(j), originalDepth);
+						  }
+						  else if (nf.node()->child(j) != BVH8::emptyNode) {
+							  BuildRecord br;
+							  br.pinfo.geomBounds = nf.node()->bounds(j);
+							  br.parent = (size_t*)&nf.node()->child(j);
+                br.depth = originalDepth;
+                //std::cout << originalDepth << std::endl;
+							  miniTreeRecordsPruned[newMiniTreeCount++] = br;
+						  }
+					  }
+				  }
+				  else if (nf != BVH8::emptyNode){
+					miniTreeRecordsPruned[newMiniTreeCount++] = miniTreeRecords[i];
+				  }
+			  }
+			  return newMiniTreeCount;
+
+		  }
+
+
+		  __forceinline void printTime(std::string s, std::chrono::high_resolution_clock::time_point& t0) const {
+
+			  std::chrono::duration<float> time = std::chrono::duration_cast<std::chrono::duration<float>>(std::chrono::high_resolution_clock::now()
+ - t0);
+			  std::cout << s << ": " << time.count() * 1000.f  << "ms" << std::endl;
+		  }
+
+		  __forceinline std::chrono::high_resolution_clock::time_point startTime() {
+			  return std::chrono::high_resolution_clock::now();
+		  }
+
+
+
+        /*! builder entry function */
+        __forceinline const ReductionTy operator() (BuildRecord& record)
+        {
+
+      // should be more sofisticated, top tree sweep is performed
+      // on an array of buildrecords, but should only be on BBoxes and indices
+      int maxMiniTrees = (numberOfTriangles/MAX_MINI_TREE_SIZE) > 0 ? (numberOfTriangles/MAX_MINI_TREE_SIZE) : 1;
+			miniTreeRecords.resize(maxMiniTrees*4);
+			miniTreeRecordsPruned.resize(maxMiniTrees*4*64);
+			miniTreeCount = 0;
+			underFilledNodes = 0;
+
+
+
+			//auto t = startTime();
+			miniTreeSelectRescursive(record);
+
+			//printTime("\nMini tree select", t);
+			//std::cout << "Mini tree count: " << miniTreeCount << std::endl;
+			miniTreeRoots.resize(miniTreeCount);
+
+		//	t = startTime();
+			buildMiniTrees(miniTreeRoots);
+		//	printTime("Build mini trees", t);
+
+		//	t = startTime();
+			if (PRUNING_THRESHOLD < 1.f && PRUNING_THRESHOLD > 0.f) {
+				int newMiniTreeCount = prune();
+        while (float(newMiniTreeCount) / float(miniTreeCount) < 10.f) {
+          PRUNING_THRESHOLD *= 0.4f;
+          newMiniTreeCount = prune();
+        }
+				miniTreeCount = newMiniTreeCount;
+			}
+		//	printTime("Pruning", t);
+			//std::cout << "Mini tree count after pruning: " << miniTreeCount << std::endl;
+		//	t = startTime();
+			if (miniTreeCount > 1)
+				buildTopTree(record);
+			else
+				record = miniTreeRecords[0];
+		//	printTime("Build top tree", t);
+			//std::cout << "Underfilled nodes: " << underFilledNodes << std::endl;
+
+			return 0;
+	    }
+
+      private:
+        Heuristic& heuristic;
+        const ReductionTy identity;
+        CreateAllocFunc& createAlloc;
+        CreateNodeFunc& createNode;
+        UpdateNodeFunc& updateNode;
+        CreateLeafFunc& createLeaf;
+        ProgressMonitor& progressMonitor;
+
+      private:
+        const PrimInfo& pinfo;
+        const size_t branchingFactor;
+        const size_t maxDepth;
+        const size_t logBlockSize;
+        const size_t minLeafSize;
+        const size_t maxLeafSize;
+        const float travCost;
+        const float intCost;
+      };
+
+    /* SAH builder that operates on an array of BuildRecords */
+    struct BVHBuilderBonsai
+    {
+      typedef range<size_t> Set;
+      typedef HeuristicBonsai<PrimRef> Heuristic;
+		typedef GeneralBuildRecord<Set, typename Heuristic::Split> BuildRecord;
+
+      /*! standard build without reduction */
+      template<typename NodeRef,
+        typename CreateAllocFunc,
+        typename CreateNodeFunc,
+        typename CreateLeafFunc,
+        typename ProgressMonitor>
+
+        static void build(NodeRef& root,
+                          CreateAllocFunc createAlloc,
+                          CreateNodeFunc createNode,
+                          CreateLeafFunc createLeaf,
+                          ProgressMonitor progressMonitor,
+                          PrimRef* prims, const PrimInfo& pinfo,
+                          const size_t branchingFactor, const size_t maxDepth, const size_t blockSize,
+                          const size_t minLeafSize, const size_t maxLeafSize,
+                          const float travCost, const float intCost, Scene* scene, const bool splitTriangles, int numOfPrimitives)
+      {
+
+        /* use dummy reduction over integers */
+        int identity = 0;
+        auto updateNode = [] (int node, int*, size_t) -> int { return 0; };
+
+
+
+        /* initiate builder */
+        build_reduce(root,
+                     createAlloc,
+                     identity,
+                     createNode,
+                     updateNode,
+                     createLeaf,
+                     progressMonitor,
+                     prims,
+                     pinfo,
+                     branchingFactor,maxDepth,blockSize,
+                     minLeafSize,maxLeafSize,travCost,intCost, scene, splitTriangles, numOfPrimitives);
+      }
+
+      /*! special builder that propagates reduction over the tree */
+      template<typename NodeRef,
+        typename CreateAllocFunc,
+        typename ReductionTy,
+        typename CreateNodeFunc,
+        typename UpdateNodeFunc,
+        typename CreateLeafFunc,
+        typename ProgressMonitor>
+
+        static ReductionTy build_reduce(NodeRef& root,
+                                        CreateAllocFunc createAlloc,
+                                        const ReductionTy& identity,
+                                        CreateNodeFunc createNode, UpdateNodeFunc updateNode, CreateLeafFunc createLeaf,
+                                        ProgressMonitor progressMonitor,
+                                        PrimRef* prims, const PrimInfo& pinfo,
+                                        const size_t branchingFactor, const size_t maxDepth, const size_t blockSize,
+                                        const size_t minLeafSize, const size_t maxLeafSize,
+                                        const float travCost, const float intCost, Scene* scene, const bool splitTriangles, int numOfPrimitives)
+      {
+        /* builder wants log2 of blockSize as input */
+        const size_t logBlockSize = __bsr(blockSize);
+        //assert((blockSize ^ (size_t(1) << logBlockSize)) == 0);
+
+        /* instantiate array binning heuristic */
+        Heuristic heuristic(prims, scene, numOfPrimitives);
+        typedef GeneralBVHBuilder<
+          BuildRecord,
+          Heuristic,
+          ReductionTy,
+          decltype(createAlloc()),
+          CreateAllocFunc,
+          CreateNodeFunc,
+          UpdateNodeFunc,
+          CreateLeafFunc,
+          ProgressMonitor> Builder;
+
+        /* instantiate builder */
+        Builder builder(heuristic,
+                        identity,
+                        createAlloc,
+                        createNode,
+                        updateNode,
+                        createLeaf,
+                        progressMonitor,
+                        pinfo,
+                        branchingFactor,maxDepth,logBlockSize,
+                        minLeafSize,maxLeafSize,travCost,intCost, numOfPrimitives, splitTriangles, scene);
+
+        /* build hierarchy */
+        //BuildRecord br(pinfo,1,(size_t*)&root,Set(0,pinfo.size()));
+        BuildRecord br(pinfo,1,(size_t*)&root,Set(0,pinfo.size()), Set(0,pinfo.size()*2));
+        return builder(br);
+      }
+    };
+
+      }
+}

kernels/xeon/builders/heuristic_bonsai.h

@@ -0,0 +1,750 @@
+// ======================================================================== //
+// Copyright 2009-2015 Intel Corporation                                    //
+//                                                                          //
+// Licensed under the Apache License, Version 2.0 (the "License");          //
+// you may not use this file except in compliance with the License.         //
+// You may obtain a copy of the License at                                  //
+//                                                                          //
+//     http://www.apache.org/licenses/LICENSE-2.0                           //
+//                                                                          //
+// Unless required by applicable law or agreed to in writing, software      //
+// distributed under the License is distributed on an "AS IS" BASIS,        //
+// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. //
+// See the License for the specific language governing permissions and      //
+// limitations under the License.                                           //
+// ======================================================================== //
+
+#pragma once
+
+#include "heuristic_binning.h"
+#include "../../algorithms/IndexSort.h"
+#include "../../common/scene.h"
+#include "presplit.h"
+#include <immintrin.h>
+
+#ifdef _WIN32
+#define ALIGNED(n) __declspec(align(n))
+#else
+#define ALIGNED(n) __attribute__((aligned(n)))
+#endif
+
+#define _MM256_TRANSPOSE4_PS(row0, row1, row2, row3) \
+do {\
+	__m256 tmp3, tmp2, tmp1, tmp0;\
+	tmp0 = _mm256_unpacklo_ps((row0), (row1));\
+	tmp2 = _mm256_unpacklo_ps((row2), (row3));\
+	tmp1 = _mm256_unpackhi_ps((row0), (row1));\
+	tmp3 = _mm256_unpackhi_ps((row2), (row3));\
+	(row0) = _mm256_movelh_ps(tmp0, tmp2);\
+	(row1) = _mm256_movehl_ps(tmp2, tmp0);\
+	(row2) = _mm256_movelh_ps(tmp1, tmp3);\
+	(row3) = _mm256_movehl_ps(tmp3, tmp1);\
+} while (0)
+
+inline __m256 combine(__m128 a, __m128 b) {
+	return _mm256_insertf128_ps(_mm256_castps128_ps256(a), b, 1);
+}
+
+inline __m256 _mm256_movelh_ps(__m256 a, __m256 b) {
+	return _mm256_shuffle_ps(a, b, _MM_SHUFFLE(1, 0, 1, 0));
+}
+
+inline __m256 _mm256_movehl_ps(__m256 a, __m256 b) {
+	return _mm256_shuffle_ps(b, a, _MM_SHUFFLE(3, 2, 3, 2));
+}
+
+inline float surfaceArea(__m256 bounds) {
+	__m128 d = _mm_add_ps(_mm256_castps256_ps128(bounds), _mm256_extractf128_ps(bounds, 1));
+	d = _mm_mul_ps(_mm_shuffle_ps(d, d, _MM_SHUFFLE(3,1,0,0)), _mm_shuffle_ps(d, d, _MM_SHUFFLE(3,2,2,1)));
+	return _mm_cvtss_f32(_mm_add_ss(_mm_add_ss(d, _mm_shuffle_ps(d, d, _MM_SHUFFLE(2,2,2,2))), _mm_shuffle_ps(d, d, _MM_SHUFFLE(1,1,1,1))));
+}
+
+namespace embree
+{
+  namespace isa
+  {
+
+
+	  struct FastSweepData {
+		  unsigned* __restrict indices[3];
+
+		  float* __restrict midX;
+		  float* __restrict midY;
+		  float* __restrict midZ;
+
+		  float* __restrict accumulatedArea;
+		  unsigned* __restrict temporaryIndices;
+		  unsigned char* __restrict leftPartition;
+		  PrimRef* __restrict gatheredRefs;
+		  unsigned offset;
+			bool isAllocated;
+
+			FastSweepData() {
+				isAllocated = false;
+			}
+
+		  FastSweepData(unsigned size, size_t maxLeafSize) {
+				fastSweepDataAllocate(size, maxLeafSize);
+		  }
+
+		  void setOffset(unsigned offset) {
+			  this->offset = offset;
+		  }
+
+			void fastSweepDataAllocate(unsigned size, size_t maxLeafSize) {
+				indices[0] = static_cast<unsigned*>(_mm_malloc(4 * size * sizeof(unsigned) + 256, 64));
+			  indices[1] = static_cast<unsigned*>(_mm_malloc(4 * size * sizeof(unsigned) + 256, 64));
+			  indices[2] = static_cast<unsigned*>(_mm_malloc(4 * size * sizeof(unsigned) + 256, 64));
+
+			  midX = static_cast<float*>(_mm_malloc(size*sizeof(float) + 256, 64));
+			  midY = static_cast<float*>(_mm_malloc(size*sizeof(float) + 256, 64));
+			  midZ = static_cast<float*>(_mm_malloc(size*sizeof(float) + 256, 64));
+
+			  accumulatedArea = static_cast<float*>(_mm_malloc(size * sizeof(float), 64)); //new float[size];
+			  temporaryIndices = new unsigned[size];
+			  leftPartition = new unsigned char[size];
+			  gatheredRefs = static_cast<PrimRef*>(_mm_malloc(maxLeafSize * 2 * sizeof(PrimRef), 64));
+				isAllocated = true;
+			}
+
+		  ~FastSweepData() {
+			  _mm_free(indices[0]);
+			  _mm_free(indices[1]);
+			  _mm_free(indices[2]);
+
+			  _mm_free(midX);
+			  _mm_free(midY);
+			  _mm_free(midZ);
+
+			  _mm_free(accumulatedArea);
+				_mm_free(gatheredRefs);
+			  delete [] temporaryIndices;
+			  delete [] leftPartition;
+		  }
+	  };
+
+	  /*! stores all information to perform some split */
+	  struct BonsaiSplit
+	  {
+		  /*! construct an invalid split by default */
+		  __forceinline BonsaiSplit()
+		  : dim(-1), pos(0.f), minimalLeafFound(false) {}
+
+		  /*! tests if this split is valid */
+		  __forceinline bool valid() const { return dim != -1; }
+
+		  /*! sets leafFound to true */
+		  __forceinline void leafFound() { minimalLeafFound = true; }
+
+		  /*! */
+		  __forceinline bool isLeaf() const { return minimalLeafFound; }
+
+		  /*! stream output */
+		  friend std::ostream& operator<<(std::ostream& cout, const BonsaiSplit& split) {
+			  return cout << "BinSplit { dim = " << split.dim << ", pos = " << split.pos << "}";
+		  }
+
+	  public:
+		  int dim;                  //!< split dimension
+		  float pos;                  //!< split plane at dim dimension
+		  int index;				//!< split index
+		  bool minimalLeafFound;	//!< true if a local SAH minimum has been found and size is less than maxLeafSize
+	  };
+
+    /*! Performs standard object binning */
+    template<typename PrimRef>
+      struct HeuristicBonsai
+      {
+		  typedef BonsaiSplit Split;
+        typedef range<size_t> Set;
+
+        static const size_t PARALLEL_THRESHOLD = 10000;
+        static const size_t PARALLEL_FIND_BLOCK_SIZE = 4096;
+        static const size_t PARALLEL_PARITION_BLOCK_SIZE = 64;
+
+				volatile int newIndex = 0;
+
+		FastSweepData** __restrict ThreadLocal;
+
+        __forceinline HeuristicBonsai ()
+          : prims(nullptr) {}
+
+        /*! remember prim array */
+        __forceinline HeuristicBonsai (PrimRef* prims, Scene* scene, int splitIndex)
+          : prims(prims), scene(scene), splitIndex(splitIndex) {}
+
+		  /*! array partitioning */
+		  void split(const Split& split, const PrimInfo& pinfo, const Set& set, PrimInfo& left, Set& lset, PrimInfo& right, Set& rset)
+		  {
+						if (likely(pinfo.size() < PARALLEL_THRESHOLD))
+					  	sequential_split(split,set,left,lset,right,rset);
+							//sequential_partition(split,set,left,lset,right,rset);
+						else
+			 				  parallel_split(split,set,left,lset,right,rset);
+		  }
+
+			void sequential_partition(const Split& split, const Set& set, PrimInfo& left, Set& lset, PrimInfo& right, Set& rset) {
+				const size_t begin = set.begin();
+			  const size_t end   = set.end();
+				CentGeomBBox3fa local_left(empty);
+				CentGeomBBox3fa local_right(empty);
+			  const float splitPos = split.pos;
+			  const unsigned int splitDim = split.dim;
+
+				int i = begin;
+				int j = end - 1;
+
+				while (prims[i].center2()[splitDim] <= splitPos) {
+					local_left.extend(prims[i].bounds());
+					++i;
+				}
+
+				while (prims[j].center2()[splitDim] > splitPos) {
+					local_right.extend(prims[j].bounds());
+					--j;
+				}
+
+				while (i < j) {
+
+					xchg(prims[i], prims[j]);
+					local_left.extend(prims[i].bounds());
+
+					local_right.extend(prims[j].bounds());
+					++i;
+					--j;
+
+					while (prims[i].center2()[splitDim] <= splitPos) {
+						local_left.extend(prims[i].bounds());
+						++i;
+					}
+
+					while (prims[j].center2()[splitDim] > splitPos) {
+						local_right.extend(prims[j].bounds());
+						--j;
+					}
+				}
+
+				size_t center = i;
+				new (&left ) PrimInfo(begin,center,local_left.geomBounds,local_left.centBounds);
+			  new (&right) PrimInfo(center,end,local_right.geomBounds,local_right.centBounds);
+			  new (&lset) range<size_t>(begin,center);
+			  new (&rset) range<size_t>(center,end);
+			}
+		  /*! array partitioning */
+		  void sequential_split(const Split& split, const Set& set, PrimInfo& left, Set& lset, PrimInfo& right, Set& rset)
+		  {
+			  /*
+			  if (!split.valid()) {
+				  deterministic_order(set);
+				  return splitFallback(set,left,lset,right,rset);
+			  }
+			  */
+
+			  const size_t begin = set.begin();
+			  const size_t end   = set.end();
+			  CentGeomBBox3fa local_left(empty);
+			  CentGeomBBox3fa local_right(empty);
+			  const float splitPos = split.pos;
+			  const unsigned int splitDim = split.dim;
+			  const unsigned int splitDimMask = (unsigned int)1 << splitDim; // FIXME: also use in unaligned and spatial binner
+
+			  size_t center = serial_partitioning(prims,begin,end,local_left,local_right,
+												  [&] (const PrimRef& ref) { return ref.center2()[splitDim] < splitPos; },
+												  [] (CentGeomBBox3fa& pinfo,const PrimRef& ref) { pinfo.extend(ref.bounds()); });
+
+			  new (&left ) PrimInfo(begin,center,local_left.geomBounds,local_left.centBounds);
+			  new (&right) PrimInfo(center,end,local_right.geomBounds,local_right.centBounds);
+			  new (&lset) range<size_t>(begin,center);
+			  new (&rset) range<size_t>(center,end);
+			  assert(area(left.geomBounds) >= 0.0f);
+			  assert(area(right.geomBounds) >= 0.0f);
+		  }
+
+		  /*! array partitioning */
+		  void parallel_split(const Split& split, const Set& set, PrimInfo& left, Set& lset, PrimInfo& right, Set& rset)
+		  {
+			  /*
+			  if (!split.valid()) {
+				  deterministic_order(set);
+				  return splitFallback(set,left,lset,right,rset);
+			  }
+			  */
+
+			  const size_t begin = set.begin();
+			  const size_t end   = set.end();
+			  left.reset();
+			  right.reset();
+			  PrimInfo init; init.reset();
+			  const float splitPos = split.pos;
+			  const unsigned int splitDim = split.dim;
+
+			  const size_t mid = parallel_in_place_partitioning<PARALLEL_PARITION_BLOCK_SIZE,PrimRef,PrimInfo>
+			  (&prims[begin],end-begin,init,left,right,
+			   [&] (const PrimRef &ref) { return ref.center2()[splitDim] < splitPos; },
+			   [] (PrimInfo &pinfo,const PrimRef &ref) { pinfo.add(ref.bounds()); },
+			   [] (PrimInfo &pinfo0,const PrimInfo &pinfo1) { pinfo0.merge(pinfo1); });
+
+			  const size_t center = begin+mid;
+			  left.begin  = begin;  left.end  = center; // FIXME: remove?
+			  right.begin = center; right.end = end;
+
+			  new (&lset) range<size_t>(begin,center);
+			  new (&rset) range<size_t>(center,end);
+		  }
+
+		  void sortIndices(FastSweepData& data, const int count, unsigned offset) {
+
+			  __m128 neg = _mm_set_ps(0.0f, -0.0f, -0.0f, -0.f); // shouldn't it be (-0, -0, -0, 0) ?
+
+
+			  int i = 0;
+			  for (; i < count - 3; i += 4) {
+
+				// crashes if prims[(i + 0 + offset)].center2() are used
+
+				//  __m128 mid0 = prims[(i + 0 + offset)].center2();
+
+				__m128 mid0 = _mm_add_ps(prims[(i + 0 + offset)].lower, prims[(i + 0 + offset)].upper);
+				  prims[(i + 0 + offset)].lower = _mm_xor_ps(prims[(i + 0 + offset)].lower, neg);
+				  //__m128 mid1 = prims[(i + 1 + data.offset)].center2();
+					__m128 mid1 = _mm_add_ps(prims[(i + 1 + offset)].lower, prims[(i + 1 + offset)].upper);
+				  prims[(i + 1 + offset)].lower = _mm_xor_ps(prims[(i + 1 + offset)].lower, neg);
+				  //__m128 mid2 = prims[(i + 2 + data.offset)].center2();
+					__m128 mid2 = _mm_add_ps(prims[(i + 2 + offset)].lower, prims[(i + 2 + offset)].upper);
+				  prims[(i + 2 + offset)].lower = _mm_xor_ps(prims[(i + 2 + offset)].lower, neg);
+				  //__m128 mid3 = prims[(i + 3 + data.offset)].center2();
+					__m128 mid3 = _mm_add_ps(prims[(i + 3 + offset)].lower, prims[(i + 3 + offset)].upper);
+				  prims[(i + 3 + offset)].lower = _mm_xor_ps(prims[(i + 3 + offset)].lower, neg);
+
+
+				  _MM_TRANSPOSE4_PS(mid0, mid1, mid2, mid3);
+				  _mm_store_ps(data.midX + i, mid0);
+				  _mm_store_ps(data.midY + i, mid1);
+				  _mm_store_ps(data.midZ + i, mid2);
+//#define STD_SORT
+
+#ifdef STD_SORT
+					data.indices[0][i+0] = data.indices[1][i+0] = data.indices[2][i+0] = i + offset + 0;
+					data.indices[0][i+1] = data.indices[1][i+1] = data.indices[2][i+1] = i + offset + 1;
+					data.indices[0][i+2] = data.indices[1][i+2] = data.indices[2][i+2] = i + offset + 2;
+					data.indices[0][i+3] = data.indices[1][i+3] = data.indices[2][i+3] = i + offset + 3;
+#endif
+			  }
+
+			  for (; i < count; i++) {
+
+				  //Vec3fa mid = prims[(i + offset)].center2();
+					__m128 mid = _mm_add_ps(prims[(i + offset)].lower, prims[(i + offset)].upper);
+				  prims[(i + offset)].lower = _mm_xor_ps(prims[(i + offset)].lower, neg);
+					ALIGNED(32) float midArray[4];
+					_mm_store_ps(midArray, mid);
+				  data.midX[i] = midArray[0];
+				  data.midY[i] = midArray[1];
+				  data.midZ[i] = midArray[2];
+#ifdef STD_SORT
+					data.indices[0][i+0] = data.indices[1][i+0] = data.indices[2][i+0] = i + offset + 0;
+#endif
+			  }
+				float* mmx = data.midX;
+				float* mmy = data.midY;
+				float* mmz = data.midZ;
+
+#ifdef STD_SORT
+				std::sort(data.indices[0], data.indices[0] + count, [mmx, offset](unsigned& a, unsigned& b) { return mmx[a-offset] < mmx[b-offset]; });
+				std::sort(data.indices[1], data.indices[1] + count, [mmy, offset](unsigned& a, unsigned& b) { return mmy[a-offset] < mmy[b-offset]; });
+				std::sort(data.indices[2], data.indices[2] + count, [mmz, offset](unsigned& a, unsigned& b) { return mmz[a-offset] < mmz[b-offset]; });
+#else
+				sortedIndicesFromFloatsX3(data.midX, data.midY, data.midZ, data.indices[0], data.indices[1], data.indices[2], count, offset);
+#endif
+				for (int i = 0; i < count-1; i++) {
+					assert(mmx[data.indices[0][i]-offset] <= mmx[data.indices[0][i+1]-offset]);
+					assert(mmy[data.indices[1][i]-offset] <= mmy[data.indices[1][i+1]-offset]);
+					assert(mmz[data.indices[2][i]-offset] <= mmz[data.indices[2][i+1]-offset]);
+				}
+			  }
+
+			__forceinline __m256 computeBounds(const float* __restrict leafBounds, const unsigned* __restrict remap, unsigned first, unsigned last) {
+				__m256 bounds0 = _mm256_load_ps(leafBounds + remap[first]*8);
+				__m256 bounds1 = bounds0;
+				__m256 bounds2 = bounds0;
+				__m256 bounds3 = bounds0;
+
+				unsigned i = first;
+
+				for (; i+7 < last; i += 8) {
+					unsigned i0 = remap[i+0];
+					unsigned i1 = remap[i+1];
+					unsigned i2 = remap[i+2];
+					unsigned i3 = remap[i+3];
+
+					__m256 a0 = _mm256_load_ps(leafBounds + i0*8);
+					__m256 a1 = _mm256_load_ps(leafBounds + i1*8);
+					__m256 a2 = _mm256_load_ps(leafBounds + i2*8);
+					__m256 a3 = _mm256_load_ps(leafBounds + i3*8);
+
+					unsigned i4 = remap[i+4];
+					unsigned i5 = remap[i+5];
+					unsigned i6 = remap[i+6];
+					unsigned i7 = remap[i+7];
+
+					__m256 a4 = _mm256_load_ps(leafBounds + i4*8);
+					__m256 a5 = _mm256_load_ps(leafBounds + i5*8);
+					__m256 a6 = _mm256_load_ps(leafBounds + i6*8);
+					__m256 a7 = _mm256_load_ps(leafBounds + i7*8);
+
+					a1 = _mm256_max_ps(a1, a0);
+					a3 = _mm256_max_ps(a3, a2);
+					a5 = _mm256_max_ps(a5, a4);
+					a7 = _mm256_max_ps(a7, a6);
+
+					bounds0 = _mm256_max_ps(bounds0, a1);
+					bounds1 = _mm256_max_ps(bounds1, a3);
+					bounds2 = _mm256_max_ps(bounds2, a5);
+					bounds3 = _mm256_max_ps(bounds3, a7);
+				}
+
+				for (; i < last; ++i)
+					bounds0 = _mm256_max_ps(bounds0, _mm256_load_ps(leafBounds + remap[i]*8));
+
+				return _mm256_max_ps(_mm256_max_ps(bounds0, bounds1), _mm256_max_ps(bounds2, bounds3));
+			}
+
+			inline BBox3fa avxBoundsToEmbreeBounds(__m256 bounds) {
+				  __m128 neg = _mm_set_ps(0.0f, -0.0f, -0.0f, -0.f);
+
+				Vec3fa lower(_mm_xor_ps(_mm256_extractf128_ps(bounds, 0), neg));
+				Vec3fa upper(_mm256_extractf128_ps(bounds, 1));
+
+				return BBox3fa(lower, upper);
+			}
+
+		  bool fastSweep(FastSweepData& data, PrimInfo& current, Split& csplit, Set& cprims, PrimInfo& left, Set& lprims, Split& lsplit, PrimInfo& right, Set& rprims, Split& rsplit, const size_t maxLeafSize, unsigned offset) {
+
+				int block_shift = 2;
+				const int mult = 1;//1000000;
+			  const float traversalCost = 0.5f;
+			  const float intersectionCost = 1.0f;
+
+			  unsigned first = cprims.begin();
+			  unsigned last = cprims.end();
+
+				 float* __restrict__ leafBounds = reinterpret_cast<float*>(prims);
+
+				/*
+				if (depth == 0) {
+					__m256 cb = computeBounds(leafBounds, data.indices[0], first, last);
+					current.geomBounds = avxBoundsToEmbreeBounds(cb);
+					float divPsa = intersectionCost / surfaceArea(cb);
+
+				}
+				*/
+			  float divPsa = intersectionCost / surfaceArea(combine(_mm_xor_ps(current.geomBounds.lower, _mm_set_ps(-0.0f, -0.0f, -0.0f, -0.f)), current.geomBounds.upper));
+
+			  float bestSah = std::numeric_limits<float>::infinity();
+
+			  unsigned bestDim = 0xffffffff;
+			  unsigned pivot = 0xffffffff;
+
+			  __m256 bestBounds;
+				//std::cout << "first: " << first << " pivot: " << pivot << " last: " << last << std::endl;
+			  if (divPsa != std::numeric_limits<float>::infinity()) {
+
+				  for (unsigned dim = 0; dim < 3; ++dim) {
+					  unsigned* __restrict sortedIndices = data.indices[dim];
+					  float* __restrict accumulatedArea = data.accumulatedArea;
+
+					  /* left to right accumulate sweep */
+					  __m256 bounds = _mm256_load_ps(leafBounds + sortedIndices[first]*8);
+
+					  {
+						  int i = (int)first;
+
+						  for(; (int)i < (int)last - 8; i+=8) {
+
+							  __m256 bb0 = _mm256_load_ps(leafBounds + sortedIndices[i+0]*8);
+							  __m256 bb1 = _mm256_load_ps(leafBounds + sortedIndices[i+1]*8);
+							  __m256 bb2 = _mm256_load_ps(leafBounds + sortedIndices[i+2]*8);
+							  __m256 bb3 = _mm256_load_ps(leafBounds + sortedIndices[i+3]*8);
+
+							  bb0 = _mm256_max_ps(bb0, bounds);
+
+							  __m256 bb4 = _mm256_load_ps(leafBounds + sortedIndices[i+4]*8);
+							  __m256 bb5 = _mm256_load_ps(leafBounds + sortedIndices[i+5]*8);
+							  __m256 bb6 = _mm256_load_ps(leafBounds + sortedIndices[i+6]*8);
+							  __m256 bb7 = _mm256_load_ps(leafBounds + sortedIndices[i+7]*8);
+
+
+							  bb1 = _mm256_max_ps(bb1, bb0);
+							  bb3 = _mm256_max_ps(bb3, bb2);
+							  bb5 = _mm256_max_ps(bb5, bb4);
+							  bb7 = _mm256_max_ps(bb7, bb6);
+							  bb3 = _mm256_max_ps(bb3, bb1);
+							  bb7 = _mm256_max_ps(bb7, bb5);
+							  bb7 = _mm256_max_ps(bb7, bb3);
+							  bb5 = _mm256_max_ps(bb5, bb3);
+							  bb4 = _mm256_max_ps(bb4, bb3);
+							  bb2 = _mm256_max_ps(bb2, bb1);
+							  bb6 = _mm256_max_ps(bb6, bb5);
+
+							  bounds = bb7;
+
+							   __m256 d0 = _mm256_add_ps(_mm256_permute2f128_ps(bb0, bb4, (0) | ((2) << 4)),
+							   _mm256_permute2f128_ps(bb0, bb4, (1) | ((3) << 4)));
+							   __m256 d1 = _mm256_add_ps(_mm256_permute2f128_ps(bb1, bb5, (0) | ((2) << 4)),
+							   _mm256_permute2f128_ps(bb1, bb5, (1) | ((3) << 4)));
+							   __m256 d2 = _mm256_add_ps(_mm256_permute2f128_ps(bb2, bb6, (0) | ((2) << 4)),
+							   _mm256_permute2f128_ps(bb2, bb6, (1) | ((3) << 4)));
+							   __m256 d3 = _mm256_add_ps(_mm256_permute2f128_ps(bb3, bb7, (0) | ((2) << 4)),
+							   _mm256_permute2f128_ps(bb3, bb7, (1) | ((3) << 4)));
+
+
+							  _MM256_TRANSPOSE4_PS(d0, d1, d2, d3);
+
+							  __m256 area = _mm256_fmadd_ps(d0, d1, _mm256_fmadd_ps(d0, d2, _mm256_mul_ps(d1, d2)));
+
+							  _mm256_storeu_ps(accumulatedArea + i, area);
+
+						  }
+							for (; i < last-1; ++i) {
+								bounds = _mm256_max_ps(bounds, _mm256_load_ps(leafBounds + sortedIndices[i]*8));
+								accumulatedArea[i] = surfaceArea(bounds);
+							}
+					  }
+
+						/* right to left sah cost computation sweep */
+						bounds = _mm256_load_ps(leafBounds + sortedIndices[last-1]*8);
+						{
+							__m256 bb[8];
+							int i = (int)last - 1;
+
+							for (; i > (int)first+7; i -= 8) {
+
+								bb[0] = _mm256_load_ps(leafBounds + sortedIndices[i-0]*8);
+								bb[1] = _mm256_load_ps(leafBounds + sortedIndices[i-1]*8);
+								bb[2] = _mm256_load_ps(leafBounds + sortedIndices[i-2]*8);
+								bb[3] = _mm256_load_ps(leafBounds + sortedIndices[i-3]*8);
+
+								bb[0] = _mm256_max_ps(bb[0], bounds);
+
+								bb[4] = _mm256_load_ps(leafBounds + sortedIndices[i-4]*8);
+								bb[5] = _mm256_load_ps(leafBounds + sortedIndices[i-5]*8);
+								bb[6] = _mm256_load_ps(leafBounds + sortedIndices[i-6]*8);
+								bb[7] = _mm256_load_ps(leafBounds + sortedIndices[i-7]*8);
+
+
+								bb[1] = _mm256_max_ps(bb[1], bb[0]);
+								bb[3] = _mm256_max_ps(bb[3], bb[2]);
+								bb[5] = _mm256_max_ps(bb[5], bb[4]);
+								bb[7] = _mm256_max_ps(bb[7], bb[6]);
+
+								bb[3] = _mm256_max_ps(bb[3], bb[1]);
+
+								__m256 lsa = _mm256_loadu_ps(accumulatedArea + (i - 8));
+
+								bb[7] = _mm256_max_ps(bb[7], bb[5]);
+
+								lsa = _mm256_shuffle_ps(lsa, lsa, _MM_SHUFFLE(0,1,2,3));
+								lsa = _mm256_permute2f128_ps(lsa, lsa, (1) | ((0) << 4));
+
+								bb[7] = _mm256_max_ps(bb[7], bb[3]);
+								bb[5] = _mm256_max_ps(bb[5], bb[3]);
+								bb[4] = _mm256_max_ps(bb[4], bb[3]);
+								bb[2] = _mm256_max_ps(bb[2], bb[1]);
+								bb[6] = _mm256_max_ps(bb[6], bb[5]);
+
+								bounds = bb[7];
+								__m256 d0 = _mm256_add_ps(_mm256_permute2f128_ps(bb[0], bb[4], (0) | ((2) << 4)),
+														  _mm256_permute2f128_ps(bb[0], bb[4], (1) | ((3) << 4)));
+								__m256 d1 = _mm256_add_ps(_mm256_permute2f128_ps(bb[1], bb[5], (0) | ((2) << 4)),
+														  _mm256_permute2f128_ps(bb[1], bb[5], (1) | ((3) << 4)));
+								__m256 d2 = _mm256_add_ps(_mm256_permute2f128_ps(bb[2], bb[6], (0) | ((2) << 4)),
+														  _mm256_permute2f128_ps(bb[2], bb[6], (1) | ((3) << 4)));
+								__m256 d3 = _mm256_add_ps(_mm256_permute2f128_ps(bb[3], bb[7], (0) | ((2) << 4)),
+														  _mm256_permute2f128_ps(bb[3], bb[7], (1) | ((3) << 4)));
+
+
+								_MM256_TRANSPOSE4_PS(d0, d1, d2, d3);
+
+								__m256 rsa = _mm256_fmadd_ps(d0, d1, _mm256_fmadd_ps(d0, d2, _mm256_mul_ps(d1, d2)));
+								__m256i leftIndices;
+								__m256i rightIndices;
+								if (last - first <= maxLeafSize*mult) {
+								int lbi[8];
+								int rbi[8];
+								for (int bi = 0; bi < 8; bi++) {
+										lbi[bi] = ((size_t((i - bi) - first)+(size_t(1)<<block_shift)-1) >> block_shift);
+										rbi[bi] = ((size_t(last - (i - bi))+(size_t(1)<<block_shift)-1) >> block_shift);
+								}
+								//std::cout << "first: " << first << " pivot: " << pivot << " last: " << last << std::endl;
+								leftIndices = _mm256_setr_epi32(lbi[0], lbi[1], lbi[2], lbi[3], lbi[4], lbi[5], lbi[6], lbi[7]);
+								rightIndices = _mm256_setr_epi32(rbi[0], rbi[1], rbi[2], rbi[3], rbi[4], rbi[5], rbi[6], rbi[7]);
+								}
+/*
+								__m256 sah = _mm256_add_ps(_mm256_mul_ps(lsa, leftIndices),
+															_mm256_mul_ps(rsa, rightIndices));
+*/
+								else {
+									__m256i indices = _mm256_sub_epi32(_mm256_set1_epi32(i), _mm256_setr_epi32(0, 1, 2, 3, 4, 5, 6, 7));
+									leftIndices = _mm256_sub_epi32(indices, _mm256_set1_epi32(first));
+									rightIndices = _mm256_sub_epi32(_mm256_set1_epi32(last), indices);
+								}
+
+/*
+								__m256 sah = _mm256_fmadd_ps(_mm256_set1_ps(intersectionCost * divPsa),
+															_mm256_fmadd_ps(lsa, _mm256_cvtepi32_ps(leftIndices),
+															_mm256_mul_ps(rsa, _mm256_cvtepi32_ps(rightIndices))), _mm256_set1_ps(traversalCost));
+*/
+															__m256 sah = _mm256_add_ps(_mm256_mul_ps(lsa, _mm256_cvtepi32_ps(leftIndices)),
+																						_mm256_mul_ps(rsa, _mm256_cvtepi32_ps(rightIndices)));
+
+				if (_mm256_movemask_ps(_mm256_cmp_ps(sah, _mm256_set1_ps(bestSah), _CMP_LT_OQ))) {
+					__m256 minSah = _mm256_min_ps(sah, _mm256_shuffle_ps(sah, sah, _MM_SHUFFLE(2,3,0,1)));
+					minSah = _mm256_min_ps(minSah, _mm256_shuffle_ps(minSah, minSah, _MM_SHUFFLE(1,0,3,2)));
+					minSah = _mm256_min_ps(minSah, _mm256_permute2f128_ps(minSah, minSah, (1) << ((0) << 4)));
+
+					unsigned mask = _mm256_movemask_ps(_mm256_cmp_ps(minSah, sah, _CMP_EQ_OQ));
+
+					ALIGNED(32) float sahArray[8];
+					_mm256_store_ps(sahArray, sah);
+
+					unsigned index = __builtin_ctz(mask);
+
+					bestBounds = bb[index];
+					bestSah = sahArray[index];
+					bestDim = dim;
+					pivot = i-(int)index;
+				}
+			}
+
+			for (; i > (int)first; --i) {
+				bounds = _mm256_max_ps(bounds, _mm256_load_ps(leafBounds + sortedIndices[i]*8));
+
+				float lsa = accumulatedArea[i - 1];
+				float rsa = surfaceArea(bounds);
+
+				if (last - first <= maxLeafSize*mult) {
+					lsa *= float((size_t(i-first)+(size_t(1)<<block_shift)-1) >> block_shift);
+					rsa *= float((size_t(last-i)+(size_t(1)<<block_shift)-1) >> block_shift);
+				}
+				else {
+					lsa *= ((float)(i - (int)first));	// Left child.
+					rsa *= ((float)((int)last - i));	// Right child.
+				}
+
+				//std::cout << lcount << ", " << rcount << std::endl;
+				float sah =  //traversalCost + intersectionCost * divPsa * // Constant cost.
+				//(lsa * ((float)(i - (int)first)) +	// Left child.
+				//rsa * ((float)((int)last - i)));	// Right child.
+				 lsa +
+				 rsa;
+
+				if (sah < bestSah) {
+					bestBounds = bounds;
+					bestSah = sah;
+					bestDim = dim;
+					pivot = i;
+				}
+			}
+					  }
+				  }
+			  }
+			  else {
+				  divPsa = 1.f;
+				  bestDim = 0;
+			  }
+				//std::cout << "first: " << first << " pivot: " << pivot << " last: " << last << std::endl;
+				float leafBlocks = (last - first);//float((size_t(last-first)+(size_t(1)<<block_shift)-1) >> block_shift);
+				if (last - first <= maxLeafSize*mult)
+					leafBlocks = float((size_t(last-first)+(size_t(1)<<block_shift)-1) >> block_shift);
+				//float blockBestSah = surfaceArea(bestBounds) * leafBlocks;
+			  //if (bestSah > leafBlocks * intersectionCost) {
+					if (traversalCost + intersectionCost * divPsa * bestSah > leafBlocks * intersectionCost) {
+				  if (last - first > maxLeafSize) {
+
+						//if (bestSah == std::numeric_limits<float>::infinity()) {
+					  	pivot = (first + last) >> 1;
+
+							bestBounds = computeBounds(leafBounds, data.indices[bestDim], pivot, last);
+						//}
+
+				  }
+
+				  else {
+					  csplit.leafFound();
+					  return true;
+				  }
+			  }
+
+			  lprims = Set(first, pivot);
+			  rprims = Set(pivot, last);
+
+
+			  partition(data, bestDim, first, last, pivot, offset);
+
+				BBox3fa finalBounds = avxBoundsToEmbreeBounds(bestBounds);
+			  right = PrimInfo(pivot, last, finalBounds, finalBounds);
+
+				bestBounds = computeBounds(leafBounds, data.indices[bestDim], first, pivot);
+				finalBounds = avxBoundsToEmbreeBounds(bestBounds);
+
+			  left = PrimInfo(first, pivot, finalBounds, finalBounds);
+
+			  return false;
+		  }
+
+		  void partition(FastSweepData& data, unsigned axis, unsigned first, unsigned last, unsigned pivot, unsigned offset) const {
+			  unsigned* __restrict ref = data.indices[axis];
+			  unsigned char* __restrict leftPartition = data.leftPartition;
+
+			  for (unsigned i = first; i < pivot; ++i)
+				  leftPartition[ref[i] - offset] = 1;
+
+			  for (unsigned i = pivot; i < last; ++i)
+				  leftPartition[ref[i] - offset] = 0;
+
+			  partitionAccordingToRef(data, (axis+1)%3, first, last, offset);
+			  partitionAccordingToRef(data, (axis+2)%3, first, last, offset);
+		  }
+
+		  void partitionAccordingToRef(FastSweepData& data, unsigned axis, unsigned first, unsigned last, unsigned offset) const {
+
+			  unsigned* __restrict temporaryIndices = data.temporaryIndices;
+			  unsigned char* __restrict leftPartition = data.leftPartition;
+
+			  unsigned* __restrict indices = data.indices[axis];
+
+			  unsigned* source = indices + first;
+			  unsigned* end = indices + last;
+
+			  unsigned* left = source;
+			  unsigned* right = temporaryIndices + first;
+
+			  while (source < end) {
+				  unsigned index = *source;
+				  unsigned l = leftPartition[index - offset];
+
+				  *left = index;
+				  *right = index;
+
+				  left += l;
+				  right += 1-l;
+
+				  ++source;
+			  }
+			  std::copy(temporaryIndices + first, right, left);
+		  }
+
+		  void gatherReferences(FastSweepData& data, const Set& current) {
+			  PrimRef* localPrims = data.gatheredRefs;
+			  for (int i = 0; i < current.size(); i++) {
+				  localPrims[i] = prims[data.indices[0][i + current.begin()]];
+			  }
+		}
+
+      private:
+        PrimRef* prims;
+				Scene* scene;
+				volatile int splitIndex; // not currently used
+      };
+  }
+}

kernels/xeon/builders/presplit.h

@@ -28,22 +28,22 @@ namespace embree
 {
   namespace isa
   {
-    __forceinline void splitTriangle(const PrimRef& prim, int dim, float pos, 
+    __forceinline void splitTriangle(const PrimRef& prim, int dim, float pos,
                                      const Vec3fa& a, const Vec3fa& b, const Vec3fa& c, PrimRef& left_o, PrimRef& right_o)
     {
       BBox3fa left = empty, right = empty;
       const Vec3fa v[3] = { a,b,c };
-      
+
       /* clip triangle to left and right box by processing all edges */
       Vec3fa v1 = v[2];
       for (size_t i=0; i<3; i++)
       {
         Vec3fa v0 = v1; v1 = v[i];
         float v0d = v0[dim], v1d = v1[dim];
-        
+
         if (v0d <= pos) left. extend(v0); // this point is on left side
         if (v0d >= pos) right.extend(v0); // this point is on right side
-        
+
         if ((v0d < pos && pos < v1d) || (v1d < pos && pos < v0d)) // the edge crosses the splitting location
         {
           assert((v1d-v0d) != 0.0f);
@@ -54,12 +54,12 @@ namespace embree
       }
       //assert(!left.empty());  // happens if split does not hit triangle
       //assert(!right.empty()); // happens if split does not hit triangle
-      
+
       /* clip against current bounds */
       BBox3fa bounds = prim.bounds();
       BBox3fa cleft (max(left .lower,bounds.lower),min(left .upper,bounds.upper));
       BBox3fa cright(max(right.lower,bounds.lower),min(right.upper,bounds.upper));
-      
+
       new (&left_o ) PrimRef(cleft, prim.geomID(), prim.primID());
       new (&right_o) PrimRef(cright,prim.geomID(), prim.primID());
     }
@@ -84,11 +84,11 @@ namespace embree
       auto order = [] (const Item& a, const Item& b) {
         return halfArea(a.sceneBounds) < halfArea(b.sceneBounds);
       };
-      
+
       while (i < N)
       {
         /* pop item with largest scene bounds */
-        std::pop_heap (&items[0],&items[i],order); 
+        std::pop_heap (&items[0],&items[i],order);
         const Item item = items[--i];
 
         /* split scene bounds into two halfes along largest dimension */
@@ -108,7 +108,7 @@ namespace embree
           items[i++] = Item(rsceneBounds,item.prim);
           std::push_heap (&items[0],&items[i],order);
         }
-        
+
         /* if split is on the right side of primitive, continue with left half of scene bounds */
         else if (center > upper)
         {
@@ -137,7 +137,7 @@ namespace embree
       for (size_t i=1; i<N; i++)
         prims_o[i-1] = items[i].prim;
     }
-    
+
     template<typename Area, typename Split>
       inline const PrimInfo presplit(const PrimInfo& pinfo, mvector<PrimRef>& prims, const Area& area, const Split& split)
     {
@@ -155,7 +155,7 @@ namespace embree
       do {
         f *= 0.9f;
         N = pinfo.size() + parallel_prefix_sum (state, size_t(0), pinfo.size(), size_t(1024), size_t(0), [&] (const range<size_t>& r, const size_t sum) -> size_t
-        { 
+        {
           size_t N=0;
           for (size_t i=r.begin(); i<r.end(); i++) {
             const float nf = ceil(f*pinfo.size()*area(prims[i])/A);
@@ -164,7 +164,7 @@ namespace embree
           }
           return N;
         },std::plus<size_t>());
-        if (iter++ == 10) // to avoid infinite loop, e.g. for pinfo.size()==1 and prims.size()==1
+        if (iter++ == 100) // to avoid infinite loop, e.g. for pinfo.size()==1 and prims.size()==1
           return pinfo;
 
       } while (pinfo.size()+N > prims.size());
@@ -209,8 +209,8 @@ namespace embree
     template<>
       inline const PrimInfo presplit<TriangleMesh>(Scene* scene, const PrimInfo& pinfo, mvector<PrimRef>& prims)
     {
-      return presplit(pinfo,prims, 
-                      [&] (const PrimRef& prim) -> float { 
+      return presplit(pinfo,prims,
+                      [&] (const PrimRef& prim) -> float {
                         const size_t geomID = prim.geomID();
                         const size_t primID = prim.primID();
                         const TriangleMesh* mesh = scene->getTriangleMesh(geomID);
@@ -240,4 +240,3 @@ namespace embree
     }
   }
 }
-

kernels/xeon/bvh8/bvh8.cpp

@@ -58,13 +58,14 @@ namespace embree
   DECLARE_BUILDER(void,Scene,size_t,BVH8Triangle4SceneBuilderSAH);
   DECLARE_BUILDER(void,Scene,size_t,BVH8Triangle8SceneBuilderSAH);
   DECLARE_BUILDER(void,Scene,size_t,BVH8TrianglePairs8SceneBuilderSAH);
+  DECLARE_BUILDER(void,Scene,size_t,BVH8Triangle4SceneBuilderBonsai);
   //DECLARE_BUILDER(void,Scene,size_t,BVH8Triangle8vSceneBuilderSAH);
 
   DECLARE_BUILDER(void,Scene,size_t,BVH8Triangle4SceneBuilderSpatialSAH);
   DECLARE_BUILDER(void,Scene,size_t,BVH8Triangle8SceneBuilderSpatialSAH);
   //DECLARE_BUILDER(void,Scene,size_t,BVH8Triangle8vSceneBuilderSpatialSAH);
 
-  void BVH8Register () 
+  void BVH8Register ()
   {
     int features = getCPUFeatures();
 
@@ -73,11 +74,11 @@ namespace embree
     SELECT_SYMBOL_AVX(features,BVH8Triangle8SceneBuilderSAH);
     //SELECT_SYMBOL_AVX(features,BVH8Triangle8vSceneBuilderSAH);
     SELECT_SYMBOL_AVX(features,BVH8TrianglePairs8SceneBuilderSAH);
-    
+
     SELECT_SYMBOL_AVX(features,BVH8Triangle4SceneBuilderSpatialSAH);
     SELECT_SYMBOL_AVX(features,BVH8Triangle8SceneBuilderSpatialSAH);
     //SELECT_SYMBOL_AVX(features,BVH8Triangle8vSceneBuilderSpatialSAH);
- 
+    SELECT_SYMBOL_AVX2(features,BVH8Triangle4SceneBuilderBonsai);
     /* select intersectors1 */
     SELECT_SYMBOL_AVX_AVX2(features,BVH8Triangle4Intersector1Moeller);
     SELECT_SYMBOL_AVX_AVX2(features,BVH8Triangle8Intersector1Moeller);
@@ -122,29 +123,29 @@ namespace embree
       numPrimitives(0), numVertices(0) {}
 
   BVH8::~BVH8 () {
-    for (size_t i=0; i<objects.size(); i++) 
+    for (size_t i=0; i<objects.size(); i++)
       delete objects[i];
   }
 
-  void BVH8::clear() 
+  void BVH8::clear()
   {
     set(BVH8::emptyNode,empty,0);
     alloc2.clear();
   }
-  
+
   void BVH8::set (NodeRef root, const BBox3fa& bounds, size_t numPrimitives)
   {
     this->root = root;
     this->bounds = bounds;
     this->numPrimitives = numPrimitives;
   }
-  
+
    void BVH8::printStatistics()
    {
      std::cout << BVH8Statistics(this).str();
      std::cout << "  "; alloc2.print_statistics();
-   }	
-  
+   }
+
   void BVH8::clearBarrier(NodeRef& node)
   {
     if (node.isBarrier())
@@ -158,7 +159,7 @@ namespace embree
 
   void BVH8::layoutLargeNodes(size_t N)
   {
-    struct NodeArea 
+    struct NodeArea
     {
       __forceinline NodeArea() {}
 
@@ -191,20 +192,20 @@ namespace embree
 
     for (size_t i=0; i<lst.size(); i++)
       lst[i].node->setBarrier();
-      
+
     root = layoutLargeNodesRecursion(root);
   }
-  
+
   BVH8::NodeRef BVH8::layoutLargeNodesRecursion(NodeRef& node)
   {
     if (node.isBarrier()) {
       node.clearBarrier();
       return node;
     }
-    else if (node.isNode()) 
+    else if (node.isNode())
     {
       Node* oldnode = node.node();
-      Node* newnode = (BVH8::Node*) alloc2.threadLocal2()->alloc0.malloc(sizeof(BVH8::Node)); // FIXME: optimize access to threadLocal2 
+      Node* newnode = (BVH8::Node*) alloc2.threadLocal2()->alloc0.malloc(sizeof(BVH8::Node)); // FIXME: optimize access to threadLocal2
       *newnode = *oldnode;
       for (size_t c=0; c<BVH8::N; c++)
         newnode->child(c) = layoutLargeNodesRecursion(oldnode->child(c));
@@ -235,7 +236,7 @@ namespace embree
     intersectors.intersector4_filter    = BVH8Triangle8Intersector4HybridMoeller;
     intersectors.intersector4_nofilter  = BVH8Triangle8Intersector4HybridMoellerNoFilter;
     intersectors.intersector8_filter    = BVH8Triangle8Intersector8HybridMoeller;
-    intersectors.intersector8_nofilter  = BVH8Triangle8Intersector8HybridMoellerNoFilter;    
+    intersectors.intersector8_nofilter  = BVH8Triangle8Intersector8HybridMoellerNoFilter;
     intersectors.intersector16_filter   = BVH8Triangle8Intersector16HybridMoeller;
     intersectors.intersector16_nofilter = BVH8Triangle8Intersector16HybridMoellerNoFilter;
     return intersectors;
@@ -270,17 +271,17 @@ namespace embree
   }
 
   Accel* BVH8::BVH8Triangle4(Scene* scene)
-  { 
+  {
     BVH8* accel = new BVH8(Triangle4::type,scene);
     Accel::Intersectors intersectors= BVH8Triangle4Intersectors(accel);
-    
+
     Builder* builder = nullptr;
     if      (scene->device->tri_builder == "default"     )          builder = BVH8Triangle4SceneBuilderSAH(accel,scene,0);
     else if (scene->device->tri_builder == "binned_sah2" )          builder = BVH8Triangle4SceneBuilderSAH(accel,scene,0);
     else if (scene->device->tri_builder == "binned_sah2_spatial" )  builder = BVH8Triangle4SceneBuilderSpatialSAH(accel,scene,0);
     else if (scene->device->tri_builder == "binned_sah2_presplit" ) builder = BVH8Triangle4SceneBuilderSAH(accel,scene,MODE_HIGH_QUALITY);
     else THROW_RUNTIME_ERROR("unknown builder "+scene->device->tri_builder+" for BVH8<Triangle4>");
-    
+
     return new AccelInstance(accel,builder,intersectors);
   }
 
@@ -300,24 +301,32 @@ namespace embree
     return new AccelInstance(accel,builder,intersectors);
   }
 
+  Accel* BVH8::BVH8Triangle4Bonsai(Scene* scene)
+	{
+		BVH8* accel = new BVH8(Triangle4::type,scene);
+		Accel::Intersectors intersectors = BVH8Triangle4Intersectors(accel);
+		Builder* builder = BVH8Triangle4SceneBuilderBonsai(accel,scene,0);
+		return new AccelInstance(accel,builder,intersectors);
+	}
+
   Accel* BVH8::BVH8Triangle8(Scene* scene)
-  { 
+  {
     BVH8* accel = new BVH8(Triangle8::type,scene);
     Accel::Intersectors intersectors= BVH8Triangle8Intersectors(accel);
-    
+
     Builder* builder = nullptr;
     if      (scene->device->tri_builder == "default"     )          builder = BVH8Triangle8SceneBuilderSAH(accel,scene,0);
     else if (scene->device->tri_builder == "binned_sah2" )          builder = BVH8Triangle8SceneBuilderSAH(accel,scene,0);
     else if (scene->device->tri_builder == "binned_sah2_spatial" )  builder = BVH8Triangle8SceneBuilderSpatialSAH(accel,scene,0);
     else if (scene->device->tri_builder == "binned_sah2_presplit" ) builder = BVH8Triangle8SceneBuilderSAH(accel,scene,MODE_HIGH_QUALITY);
     else THROW_RUNTIME_ERROR("unknown builder "+scene->device->tri_builder+" for BVH8<Triangle8>");
-    
+
     return new AccelInstance(accel,builder,intersectors);
   }
 
 
   /*Accel* BVH8::BVH8Triangle8v(Scene* scene)
-  { 
+  {
     BVH8* accel = new BVH8(Triangle8v::type,scene);
     Accel::Intersectors intersectors= BVH8Triangle8vIntersectors(accel);
     Builder* builder = nullptr;
@@ -325,20 +334,20 @@ namespace embree
     else if (scene->device->tri_builder == "binned_sah2" )          builder = BVH8Triangle8vSceneBuilderSAH(accel,scene,0);
     else if (scene->device->tri_builder == "binned_sah2_spatial" )  builder = BVH8Triangle8vSceneBuilderSpatialSAH(accel,scene,0);
     else if (scene->device->tri_builder == "binned_sah2_presplit" ) builder = BVH8Triangle8vSceneBuilderSAH(accel,scene,MODE_HIGH_QUALITY);
-    else 
-      THROW_RUNTIME_ERROR("unknown builder "+scene->device->tri_builder+" for BVH8<Triangle8v>");    
+    else
+      THROW_RUNTIME_ERROR("unknown builder "+scene->device->tri_builder+" for BVH8<Triangle8v>");
     return new AccelInstance(accel,builder,intersectors);
     }*/
 
   Accel* BVH8::BVH8TrianglePairs8(Scene* scene)
-  { 
+  {
     BVH8* accel = new BVH8(TrianglePairs8::type,scene);
     Accel::Intersectors intersectors = BVH8TrianglePairs8Intersectors(accel);
-    
+
     Builder* builder = nullptr;
     if      (scene->device->tri_builder == "default"     ) builder = BVH8TrianglePairs8SceneBuilderSAH(accel,scene,0);
     else THROW_RUNTIME_ERROR("unknown builder "+scene->device->tri_builder+" for BVH8<Triangle8>");
-    
+
     return new AccelInstance(accel,builder,intersectors);
   }
 
@@ -375,4 +384,3 @@ namespace embree
     }*/
 
 }
-

kernels/xeon/bvh8/bvh8.h

@@ -31,7 +31,7 @@ namespace embree
   {
     ALIGNED_CLASS;
   public:
-    
+
     /*! forward declaration of node type */
     struct Node;
 
@@ -44,20 +44,20 @@ namespace embree
     static const size_t alignment = 5;
 
     /*! Masks the bits that store the number of items per leaf. */
-    static const size_t align_mask = (1 << alignment)-1;  
-    static const size_t items_mask = (1 << (alignment-1))-1;  
+    static const size_t align_mask = (1 << alignment)-1;
+    static const size_t items_mask = (1 << (alignment-1))-1;
 
     /*! Empty node */
     static const size_t emptyNode = 1;
 
     /*! Invalid node, used as marker in traversal */
     static const size_t invalidNode = (((size_t)-1) & (~items_mask)) | 1;
-      
+
     /*! Maximal depth of the BVH. */
     static const size_t maxBuildDepth = 32;
     static const size_t maxBuildDepthLeaf = maxBuildDepth+16;
     static const size_t maxDepth = maxBuildDepthLeaf+maxBuildDepthLeaf+maxBuildDepth;
-    
+
     /*! Maximal number of primitive blocks in a leaf. */
     static const size_t maxLeafBlocks = 6; //items_mask-1;
 
@@ -94,7 +94,7 @@ namespace embree
 
       /*! Sets the barrier bit. */
       __forceinline void setBarrier() { ptr |= (size_t)(1 << (alignment-1)); }
-      
+
       /*! Clears the barrier bit. */
       __forceinline void clearBarrier() { ptr &= ~(size_t)(1 << (alignment-1)); }
 
@@ -103,14 +103,14 @@ namespace embree
 
       /*! checks if this is a leaf */
       __forceinline int isLeaf() const { return (ptr & (size_t)align_mask) != 0; }
-      
+
       /*! checks if this is a node */
       __forceinline int isNode() const { return (ptr & (size_t)align_mask) == 0; }
-      
+
       /*! returns node pointer */
       __forceinline       Node* node()       { assert(isNode()); return (      Node*)ptr; }
       __forceinline const Node* node() const { assert(isNode()); return (const Node*)ptr; }
-      
+
       /*! returns leaf pointer */
       __forceinline char* leaf(size_t& num) const {
         assert(isLeaf());
@@ -129,13 +129,13 @@ namespace embree
     {
       /*! Clears the node. */
       __forceinline void clear() {
-        lower_x = lower_y = lower_z = pos_inf; 
+        lower_x = lower_y = lower_z = pos_inf;
         upper_x = upper_y = upper_z = neg_inf;
 	for (size_t i=0; i<N; i++) children[i] = emptyNode;
       }
 
       /*! Sets bounding box of child. */
-      __forceinline void set(size_t i, const BBox3fa& bounds) 
+      __forceinline void set(size_t i, const BBox3fa& bounds)
       {
         assert(i < N);
         lower_x[i] = bounds.lower.x; lower_y[i] = bounds.lower.y; lower_z[i] = bounds.lower.z;
@@ -156,7 +156,7 @@ namespace embree
       }
 
       /*! Returns bounds of specified child. */
-      __forceinline BBox3fa bounds(size_t i) const 
+      __forceinline BBox3fa bounds(size_t i) const
       {
         assert(i < N);
         const Vec3fa lower(lower_x[i],lower_y[i],lower_z[i]);
@@ -219,7 +219,7 @@ namespace embree
       }
     }
 #endif
-    
+
   public:
 
     /*! BVH8 default constructor. */
@@ -232,7 +232,7 @@ namespace embree
     static Accel* BVH8Triangle4(Scene* scene);
     static Accel* BVH8Triangle4ObjectSplit(Scene* scene);
     static Accel* BVH8Triangle4SpatialSplit(Scene* scene);
-
+    static Accel* BVH8Triangle4Bonsai(Scene* scene);
     static Accel* BVH8Triangle8(Scene* scene);
     //static Accel* BVH8Triangle8v(Scene* scene);
     static Accel* BVH8TrianglePairs8(Scene* scene);
@@ -260,19 +260,19 @@ namespace embree
 #if defined (__AVX__)
 
     /*! Encodes a node */
-    __forceinline NodeRef encodeNode(Node* node) { 
+    __forceinline NodeRef encodeNode(Node* node) {
       return NodeRef((size_t) node);
     }
-    
+
     /*! Encodes a leaf */
     __forceinline NodeRef encodeLeaf(void* tri, size_t num) {
-      assert(!((size_t)tri & align_mask)); 
+      assert(!((size_t)tri & align_mask));
       return NodeRef((size_t)tri | (1+min(num,(size_t)maxLeafBlocks)));
     }
 #endif
 
   public:
-    
+
     /*! calculates the amount of bytes allocated */
     size_t bytesAllocated() {
       return alloc2.getAllocatedBytes();
@@ -305,7 +305,7 @@ namespace embree
     for (size_t i=0;i<BVH8::N;i++) o << n.children[i] << " ";
     o << std::endl;
     return o;
-  } 
+  }
 #endif
 
 }

kernels/xeon/bvh8/bvh8_builder_bonsai.avx.cpp

@@ -0,0 +1,25 @@
+// ======================================================================== //
+// Copyright 2009-2015 Intel Corporation                                    //
+//                                                                          //
+// Licensed under the Apache License, Version 2.0 (the "License");          //
+// you may not use this file except in compliance with the License.         //
+// You may obtain a copy of the License at                                  //
+//                                                                          //
+//     http://www.apache.org/licenses/LICENSE-2.0                           //
+//                                                                          //
+// Unless required by applicable law or agreed to in writing, software      //
+// distributed under the License is distributed on an "AS IS" BASIS,        //
+// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. //
+// See the License for the specific language governing permissions and      //
+// limitations under the License.                                           //
+// ======================================================================== //
+
+// We cannot compile the same file containing lambda functions for two
+// ISAs, as a lambda name mangling bug of ICC under Windows causes
+// symbols to conflict.
+
+
+//created by Per Ganestam
+
+#include "bvh8_builder_bonsai.cpp"
+

kernels/xeon/bvh8/bvh8_builder_bonsai.cpp

@@ -0,0 +1,285 @@
+// ======================================================================== //
+// Copyright 2009-2015 Intel Corporation                                    //
+//                                                                          //
+// Licensed under the Apache License, Version 2.0 (the "License");          //
+// you may not use this file except in compliance with the License.         //
+// You may obtain a copy of the License at                                  //
+//                                                                          //
+//     http://www.apache.org/licenses/LICENSE-2.0                           //
+//                                                                          //
+// Unless required by applicable law or agreed to in writing, software      //
+// distributed under the License is distributed on an "AS IS" BASIS,        //
+// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. //
+// See the License for the specific language governing permissions and      //
+// limitations under the License.                                           //
+// ======================================================================== //
+
+#include "bvh8.h"
+#include "bvh8_statistics.h"
+//#include "common/profile.h"
+
+#include "../builders/primrefgen.h"
+#include "../builders/presplit.h"
+#include "../builders/bvh_builder_bonsai.h"
+
+#include "../geometry/triangle4.h"
+#include "../geometry/triangle8.h"
+
+#define PROFILE 0
+
+namespace embree
+{
+  namespace isa
+  {
+    typedef FastAllocator::ThreadLocal2 Allocator;
+
+    struct CreateBVH8Alloc
+    {
+      __forceinline CreateBVH8Alloc (BVH8* bvh) : bvh(bvh) {}
+      __forceinline Allocator* operator() () const { return bvh->alloc2.threadLocal2();  }
+
+      BVH8* bvh;
+    };
+
+    struct CreateBVH8Node
+    {
+      __forceinline CreateBVH8Node (BVH8* bvh) : bvh(bvh) {}
+
+      __forceinline int operator() (const isa::BVHBuilderBonsai::BuildRecord& current, BVHBuilderBonsai::BuildRecord* children, const size_t N, Allocator* alloc)
+      {
+        BVH8::Node* node = nullptr;
+        node = (BVH8::Node*) alloc->alloc0.malloc(sizeof(BVH8::Node), 1 << BVH8::alignment);
+        node->clear();
+        for (size_t i=0; i<N; i++) {
+          node->set(i,children[i].pinfo.geomBounds);
+          children[i].parent = (size_t*) &node->child(i);
+        }
+        *current.parent = bvh->encodeNode(node);
+
+		return 0;
+      }
+
+		__forceinline int operator() (const isa::BVHBuilderBonsai::BuildRecord& current, BVHBuilderBonsai::BuildRecord* children, const size_t N, Allocator* alloc, int& underFilledNodes)
+		{
+			BVH8::Node* node = nullptr;
+			node = (BVH8::Node*) alloc->alloc0.malloc(sizeof(BVH8::Node), 1 << BVH8::alignment);
+			node->clear();
+			for (size_t i=0; i<N; i++) {
+				node->set(i,children[i].pinfo.geomBounds);
+				node->child(i) = BVH8::NodeRef(*children[i].parent);
+			}
+			*current.parent = bvh->encodeNode(node);
+			return 0;
+		}
+
+      BVH8* bvh;
+    };
+
+    template<typename Primitive>
+    struct CreateBVH8Leaf
+    {
+      __forceinline CreateBVH8Leaf (BVH8* bvh, PrimRef* prims) : bvh(bvh), prims(prims) {}
+
+      __forceinline int operator() (const BVHBuilderBonsai::BuildRecord& current, Allocator* alloc)
+      {
+        size_t items = Primitive::blocks(current.prims.size());
+        size_t start = current.prims.begin();
+        Primitive* accel = (Primitive*) alloc->alloc1.malloc(items*sizeof(Primitive), 1 << BVH8::alignment);
+        BVH8::NodeRef node = bvh->encodeLeaf((char*)accel,items);
+        for (size_t i=0; i<items; i++) {
+          accel[i].fill(prims,start,current.prims.end(),bvh->scene,false);
+        }
+        *current.parent = node;
+	return 1;
+      }
+
+		__forceinline int operator() (const BVHBuilderBonsai::BuildRecord& current, PrimRef* gatheredRefs, Allocator* alloc)
+		{
+			size_t items = Primitive::blocks(current.prims.size());
+			size_t start = 0; //< gathered refs start at 0 instead of current.prims.begin();
+			Primitive* accel = (Primitive*) alloc->alloc1.malloc(items*sizeof(Primitive), 1 << BVH8::alignment);
+			BVH8::NodeRef node = bvh->encodeLeaf((char*)accel,items);
+			for (size_t i=0; i<items; i++) {
+				accel[i].fill(gatheredRefs,start,current.prims.size(),bvh->scene,false); //< last ref is prims.size() instead of prims.end()
+			}
+			*current.parent = node;
+			return 1;
+		}
+
+      BVH8* bvh;
+      PrimRef* prims;
+    };
+
+    /************************************************************************************/
+    /************************************************************************************/
+    /************************************************************************************/
+    /************************************************************************************/
+
+    template<typename Mesh, typename Primitive>
+    struct BVH8BuilderBonsai : public Builder
+    {
+      BVH8* bvh;
+      Scene* scene;
+      Mesh* mesh;
+      mvector<PrimRef> prims;
+      const size_t sahBlockSize;
+      const float intCost;
+      const size_t minLeafSize;
+      const size_t maxLeafSize;
+      const float presplitFactor;
+
+      BVH8BuilderBonsai (BVH8* bvh, Scene* scene, const size_t leafBlockSize, const size_t sahBlockSize, const float intCost, const size_t minLeafSize, const size_t maxLeafSize, const size_t mode)
+        : bvh(bvh), scene(scene), mesh(nullptr), prims(scene->device), sahBlockSize(sahBlockSize), intCost(intCost), minLeafSize(minLeafSize), maxLeafSize(min(maxLeafSize,leafBlockSize*BVH8::maxLeafBlocks)),
+		presplitFactor((mode & MODE_HIGH_QUALITY) ? 1.5f : 1.0f) {}
+
+      BVH8BuilderBonsai (BVH8* bvh, Mesh* mesh, const size_t leafBlockSize, const size_t sahBlockSize, const float intCost, const size_t minLeafSize, const size_t maxLeafSize, const size_t mode)
+        : bvh(bvh), scene(nullptr), mesh(mesh), prims(scene->device), sahBlockSize(sahBlockSize), intCost(intCost), minLeafSize(minLeafSize), maxLeafSize(min(maxLeafSize,leafBlockSize*BVH8::maxLeafBlocks)),
+		presplitFactor((mode & MODE_HIGH_QUALITY) ? 1.5f : 1.0f) {}
+
+      void build(size_t, size_t)
+      {
+	/* skip build for empty scene */
+	const size_t numPrimitives = mesh ? mesh->size() : scene->getNumPrimitives<Mesh,1>();
+        if (numPrimitives == 0) {
+          prims.resize(numPrimitives);
+          bvh->set(BVH8::emptyNode,empty,0);
+          return;
+        }
+        const size_t numSplitPrimitives = max(numPrimitives,size_t(presplitFactor*numPrimitives));
+
+        /* verbose mode */
+        if (bvh->device->verbosity(1) && mesh == nullptr)
+	  std::cout << "building BVH8<" << bvh->primTy.name << "> with " << TOSTRING(isa) "::BVH8BuilderBonsai " << (presplitFactor != 1.0f ? "presplit" : "") << " ... " << std::flush;
+
+    bool splitTriangles = false; // not currently used
+
+	double t0 = 0.0f, dt = 0.0f;
+#if PROFILE
+	profile(2,20,numPrimitives,[&] (ProfileTimer& timer)
+        {
+#endif
+
+          if ((bvh->device->benchmark || bvh->device->verbosity(1)) && mesh == nullptr) t0 = getSeconds();
+
+          auto progress = [&] (size_t dn) { bvh->scene->progressMonitor(dn); };
+          auto virtualprogress = BuildProgressMonitorFromClosure(progress);
+
+	    bvh->alloc2.init_estimate(numSplitPrimitives*sizeof(PrimRef));
+	    if (splitTriangles)
+        prims.resize(numSplitPrimitives*4);
+      else
+        prims.resize(numSplitPrimitives);
+
+      PrimInfo pinfo = mesh ? createPrimRefArray<Mesh>(mesh,prims,virtualprogress) : createPrimRefArray<Mesh,1>(scene,prims,virtualprogress);
+            if (presplitFactor > 1.0f)
+              pinfo = presplit<Mesh>(scene, pinfo, prims);
+	    BVH8::NodeRef root;
+            BVHBuilderBonsai::build<BVH8::NodeRef>
+              (root,CreateBVH8Alloc(bvh),CreateBVH8Node(bvh),CreateBVH8Leaf<Primitive>(bvh,prims.data()), progress,
+               prims.data(),pinfo,BVH8::N,BVH8::maxBuildDepthLeaf,sahBlockSize,minLeafSize,maxLeafSize,BVH8::travCost,intCost, bvh->scene, splitTriangles, pinfo.size());
+
+            bvh->set(root,pinfo.geomBounds,pinfo.size());
+            //bvh->layoutLargeNodes(numSplitPrimitives*0.005f);
+
+	    if ((bvh->device->benchmark || bvh->device->verbosity(1)) && mesh == nullptr) dt = getSeconds()-t0;
+
+#if PROFILE
+           dt = timer.avg();
+        });
+#endif
+
+//  exit(0);
+	/* clear temporary data for static geometry */
+	bool staticGeom = mesh ? mesh->isStatic() : scene->isStatic();
+	if (staticGeom) prims.clear();
+	bvh->alloc2.cleanup();
+
+	/* verbose mode */
+	if (bvh->device->verbosity(1) && mesh == nullptr)
+	  std::cout << "[DONE] " << 1000.0f*dt << "ms (" << numPrimitives/dt*1E-6 << " Mtris/s)" << std::endl;
+	if (bvh->device->verbosity(2) && mesh == nullptr)
+	  bvh->printStatistics();
+
+        /* benchmark mode */
+          if (bvh->device->benchmark) {
+          BVH8Statistics stat(bvh);
+          std::cout << "BENCHMARK_BUILD " << dt << " " << double(numPrimitives)/dt << " " << stat.sah() << " " << stat.bytesUsed() << std::endl;
+        }
+      }
+
+      void clear() {
+        prims.clear();
+      }
+    };
+
+    /* entry functions for the scene builder */
+    Builder* BVH8Triangle4SceneBuilderBonsai  (void* bvh, Scene* scene, size_t mode) { return new BVH8BuilderBonsai<TriangleMesh,Triangle4>((BVH8*)bvh,scene,4,4,1.0f,4,inf,mode); }
+    Builder* BVH8Triangle8SceneBuilderBonsai  (void* bvh, Scene* scene, size_t mode) { return new BVH8BuilderBonsai<TriangleMesh,Triangle8>((BVH8*)bvh,scene,8,4,1.0f,8,inf,mode); }
+
+    /************************************************************************************/
+    /************************************************************************************/
+    /************************************************************************************/
+    /************************************************************************************/
+#if 0
+    struct CreateListBVH8Node // FIXME: merge with above class
+    {
+      __forceinline CreateListBVH8Node (BVH8* bvh) : bvh(bvh) {}
+
+      __forceinline BVH8::Node* operator() (const isa::BVHBuilderBinnedSpatialSAH::BuildRecord& current, BVHBuilderBinnedSpatialSAH::BuildRecord* children, const size_t N, Allocator* alloc)
+      {
+        BVH8::Node* node = (BVH8::Node*) alloc->alloc0.malloc(sizeof(BVH8::Node), 1 << BVH8::alignment); node->clear();
+        for (size_t i=0; i<N; i++) {
+          node->set(i,children[i].pinfo.geomBounds);
+          children[i].parent = (size_t*) &node->child(i);
+        }
+        *current.parent = bvh->encodeNode(node);
+	return node;
+      }
+
+      BVH8* bvh;
+    };
+
+
+    template<typename Primitive>
+    struct CreateBVH8ListLeaf
+    {
+      __forceinline CreateBVH8ListLeaf (BVH8* bvh) : bvh(bvh) {}
+
+      __forceinline size_t operator() (BVHBuilderBinnedSpatialSAH::BuildRecord& current, Allocator* alloc) // FIXME: why are prims passed here but not for createNode
+      {
+        size_t n = current.pinfo.size();
+        size_t N = Primitive::blocks(n);
+        Primitive* leaf = (Primitive*) alloc->alloc1.malloc(N*sizeof(Primitive), 1 << BVH8::alignment);
+        BVH8::NodeRef node = bvh->encodeLeaf((char*)leaf,N);
+
+        /*PrimRefList::block_iterator_unsafe iter1(current.prims);
+        while (iter1) {
+          const int i = iter1->lower.a;
+          iter1->lower.a = geomIDprimID[i].first;
+          iter1->upper.a = geomIDprimID[i].second;
+          iter1++;
+          }*/
+        PrimRefList::block_iterator_unsafe iter1(current.prims);
+        while (iter1) {
+          iter1->lower.a &= 0x00FFFFFF; // FIXME: hack
+          iter1++;
+        }
+
+        /* insert all triangles */
+        PrimRefList::block_iterator_unsafe iter(current.prims);
+        for (size_t i=0; i<N; i++) leaf[i].fill(iter,bvh->scene,false);
+        assert(!iter);
+
+        /* free all primitive blocks */
+        while (PrimRefList::item* block = current.prims.take())
+          delete block;
+
+        *current.parent = node;
+	return n;
+      }
+
+      BVH8* bvh;
+    };
+#endif
+  }
+}

kernels/xeon/bvh8/bvh8_builder_sah.cpp

@@ -26,7 +26,7 @@
 #include "../geometry/triangle8v.h"
 #include "../geometry/trianglepairs8.h"
 
-#define PROFILE 0
+#define PROFILE 1
 
 namespace embree
 {
@@ -45,13 +45,13 @@ namespace embree
     struct CreateBVH8Node
     {
       __forceinline CreateBVH8Node (BVH8* bvh) : bvh(bvh) {}
-      
-      __forceinline int operator() (const isa::BVHBuilderBinnedSAH::BuildRecord& current, BVHBuilderBinnedSAH::BuildRecord* children, const size_t N, Allocator* alloc) 
+
+      __forceinline int operator() (const isa::BVHBuilderBinnedSAH::BuildRecord& current, BVHBuilderBinnedSAH::BuildRecord* children, const size_t N, Allocator* alloc)
       {
         BVH8::Node* node = nullptr;
         //if (current.pinfo.size() > 4096) node = (BVH8::Node*)   bvh->alloc2.malloc(sizeof(BVH8::Node),sizeof(BVH8::Node));
         //else
-        node = (BVH8::Node*) alloc->alloc0.malloc(sizeof(BVH8::Node), 1 << BVH8::alignment); 
+        node = (BVH8::Node*) alloc->alloc0.malloc(sizeof(BVH8::Node), 1 << BVH8::alignment);
         node->clear();
         for (size_t i=0; i<N; i++) {
           node->set(i,children[i].pinfo.geomBounds);
@@ -68,7 +68,7 @@ namespace embree
     struct CreateBVH8Leaf
     {
       __forceinline CreateBVH8Leaf (BVH8* bvh, PrimRef* prims) : bvh(bvh), prims(prims) {}
-      
+
       __forceinline int operator() (const BVHBuilderBinnedSAH::BuildRecord& current, Allocator* alloc)
       {
         size_t items = Primitive::blocks(current.prims.size());
@@ -85,8 +85,8 @@ namespace embree
       BVH8* bvh;
       PrimRef* prims;
     };
-    
-    /************************************************************************************/ 
+
+    /************************************************************************************/
     /************************************************************************************/
     /************************************************************************************/
     /************************************************************************************/
@@ -106,13 +106,13 @@ namespace embree
 
       BVH8BuilderSAH (BVH8* bvh, Scene* scene, const size_t leafBlockSize, const size_t sahBlockSize, const float intCost, const size_t minLeafSize, const size_t maxLeafSize, const size_t mode)
         : bvh(bvh), scene(scene), mesh(nullptr), prims(scene->device), sahBlockSize(sahBlockSize), intCost(intCost), minLeafSize(minLeafSize), maxLeafSize(min(maxLeafSize,leafBlockSize*BVH8::maxLeafBlocks)),
-          presplitFactor((mode & MODE_HIGH_QUALITY) ? 1.5f : 1.0f) {}
+          presplitFactor((mode & MODE_HIGH_QUALITY) ? 2.5f : 1.0f) {}
 
       BVH8BuilderSAH (BVH8* bvh, Mesh* mesh, const size_t leafBlockSize, const size_t sahBlockSize, const float intCost, const size_t minLeafSize, const size_t maxLeafSize, const size_t mode)
         : bvh(bvh), scene(nullptr), mesh(mesh), prims(bvh->device), sahBlockSize(sahBlockSize), intCost(intCost), minLeafSize(minLeafSize), maxLeafSize(min(maxLeafSize,leafBlockSize*BVH8::maxLeafBlocks)),
-          presplitFactor((mode & MODE_HIGH_QUALITY) ? 1.5f : 1.0f) {}
+          presplitFactor((mode & MODE_HIGH_QUALITY) ? 2.5f : 1.0f) {}
 
-      void build(size_t, size_t) 
+      void build(size_t, size_t)
       {
 	/* skip build for empty scene */
 	const size_t numPrimitives = mesh ? mesh->size() : scene->getNumPrimitives<Mesh,1>();
@@ -122,7 +122,7 @@ namespace embree
           return;
         }
         const size_t numSplitPrimitives = max(numPrimitives,size_t(presplitFactor*numPrimitives));
-      
+
         /* verbose mode */
         if (bvh->device->verbosity(1) && mesh == nullptr)
 	  std::cout << "building BVH8<" << bvh->primTy.name << "> with " << TOSTRING(isa) "::BVH8BuilderSAH " << (presplitFactor != 1.0f ? "presplit" : "") << " ... " << std::flush;
@@ -132,18 +132,18 @@ namespace embree
 	profile(2,20,numPrimitives,[&] (ProfileTimer& timer)
         {
 #endif
-	    
+
           if ((bvh->device->benchmark || bvh->device->verbosity(1)) && mesh == nullptr) t0 = getSeconds();
-          
+
           auto progress = [&] (size_t dn) { bvh->scene->progressMonitor(dn); };
           auto virtualprogress = BuildProgressMonitorFromClosure(progress);
-          
+
 	    bvh->alloc2.init_estimate(numSplitPrimitives*sizeof(PrimRef));
 	    prims.resize(numSplitPrimitives);
 	    PrimInfo pinfo = mesh ? createPrimRefArray<Mesh>(mesh,prims,virtualprogress) : createPrimRefArray<Mesh,1>(scene,prims,virtualprogress);
             if (presplitFactor > 1.0f)
               pinfo = presplit<Mesh>(scene, pinfo, prims);
-	    BVH8::NodeRef root; 
+	    BVH8::NodeRef root;
             BVHBuilderBinnedSAH::build<BVH8::NodeRef>
               (root,CreateBVH8Alloc(bvh),CreateBVH8Node(bvh),CreateBVH8Leaf<Primitive>(bvh,prims.data()), progress,
                prims.data(),pinfo,BVH8::N,BVH8::maxBuildDepthLeaf,sahBlockSize,minLeafSize,maxLeafSize,BVH8::travCost,intCost);
@@ -155,8 +155,8 @@ namespace embree
 
 #if PROFILE
            dt = timer.avg();
-        }); 
-#endif	
+        });
+#endif
 
 	/* clear temporary data for static geometry */
 	bool staticGeom = mesh ? mesh->isStatic() : scene->isStatic();
@@ -193,7 +193,7 @@ namespace embree
     Builder* BVH8TrianglePairs8SceneBuilderSAH  (void* bvh, Scene* scene, size_t mode) { return new BVH8BuilderSAH<TriangleMesh,Triangle8>((BVH8*)bvh,scene,8,4,1.0f,8,inf,mode); }
 
 
-    /************************************************************************************/ 
+    /************************************************************************************/
     /************************************************************************************/
     /************************************************************************************/
     /************************************************************************************/
@@ -201,8 +201,8 @@ namespace embree
     struct CreateListBVH8Node // FIXME: merge with above class
     {
       __forceinline CreateListBVH8Node (BVH8* bvh) : bvh(bvh) {}
-      
-      __forceinline BVH8::Node* operator() (const isa::BVHBuilderBinnedSpatialSAH::BuildRecord& current, BVHBuilderBinnedSpatialSAH::BuildRecord* children, const size_t N, Allocator* alloc) 
+
+      __forceinline BVH8::Node* operator() (const isa::BVHBuilderBinnedSpatialSAH::BuildRecord& current, BVHBuilderBinnedSpatialSAH::BuildRecord* children, const size_t N, Allocator* alloc)
       {
         BVH8::Node* node = (BVH8::Node*) alloc->alloc0.malloc(sizeof(BVH8::Node), 1 << BVH8::alignment); node->clear();
         for (size_t i=0; i<N; i++) {
@@ -220,7 +220,7 @@ namespace embree
     struct CreateBVH8ListLeaf
     {
       __forceinline CreateBVH8ListLeaf (BVH8* bvh) : bvh(bvh) {}
-      
+
       __forceinline size_t operator() (BVHBuilderBinnedSpatialSAH::BuildRecord& current, Allocator* alloc) // FIXME: why are prims passed here but not for createNode
       {
         size_t n = current.pinfo.size();
@@ -245,7 +245,7 @@ namespace embree
         PrimRefList::block_iterator_unsafe iter(current.prims);
         for (size_t i=0; i<N; i++) leaf[i].fill(iter,bvh->scene,false);
         assert(!iter);
-        
+
         /* free all primitive blocks */
         while (PrimRefList::item* block = current.prims.take())
           delete block;
@@ -277,7 +277,7 @@ namespace embree
         : bvh(bvh), scene(nullptr), mesh(mesh), sahBlockSize(sahBlockSize), intCost(intCost), minLeafSize(minLeafSize), maxLeafSize(min(maxLeafSize,leafBlockSize*BVH8::maxLeafBlocks)),
           presplitFactor((mode & MODE_HIGH_QUALITY) ? 1.5f : 1.0f) {}
 
-      void build(size_t, size_t) 
+      void build(size_t, size_t)
       {
 	/* skip build for empty scene */
 	const size_t numPrimitives = mesh ? mesh->size() : scene->getNumPrimitives<Mesh,1>();
@@ -287,20 +287,20 @@ namespace embree
           return;
         }
         const size_t numSplitPrimitives = max(numPrimitives,size_t(presplitFactor*numPrimitives));
-      
+
         /* reduction function */
 	auto rotate = [&] (BVH8::Node* node, const size_t* counts, const size_t N) -> size_t
 	{
           size_t n = 0;
 #if ROTATE_TREE
 	  assert(N <= BVH8::N);
-          for (size_t i=0; i<N; i++) 
+          for (size_t i=0; i<N; i++)
             n += counts[i];
           if (n >= 4096) {
             for (size_t i=0; i<N; i++) {
               if (counts[i] < 4096) {
-                for (int j=0; j<ROTATE_TREE; j++) 
-                  BVH8Rotate::rotate(bvh,node->child(i)); 
+                for (int j=0; j<ROTATE_TREE; j++)
+                  BVH8Rotate::rotate(bvh,node->child(i));
                 node->child(i).setBarrier();
               }
             }
@@ -318,9 +318,9 @@ namespace embree
 	profile(2,20,numPrimitives,[&] (ProfileTimer& timer)
         {
 #endif
-	    
+
           if ((bvh->device->benchmark || bvh->device->verbosity(1)) && mesh == nullptr) t0 = getSeconds();
-	    
+
             auto progress = [&] (size_t dn) { bvh->scene->progressMonitor(dn); };
             auto virtualprogress = BuildProgressMonitorFromClosure(progress);
 
@@ -340,7 +340,7 @@ namespace embree
             {
               double A = 0.0f;
               while (PrimRefList::item* block = iter.next()) {
-                for (size_t i=0; i<block->size(); i++) 
+                for (size_t i=0; i<block->size(); i++)
                   A += area(block->at(i).bounds());
                 //A += heuristic(block->at(i));
               }
@@ -372,7 +372,7 @@ namespace embree
 	      (root,CreateBVH8Alloc(bvh),size_t(0),CreateListBVH8Node(bvh),rotate,CreateBVH8ListLeaf<Primitive>(bvh),
                [&] (const PrimRef& prim, int dim, float pos, PrimRef& left_o, PrimRef& right_o)
                {
-                TriangleMesh* mesh = (TriangleMesh*) scene->get(prim.geomID() & 0x00FFFFFF); 
+                TriangleMesh* mesh = (TriangleMesh*) scene->get(prim.geomID() & 0x00FFFFFF);
                 TriangleMesh::Triangle tri = mesh->triangle(prim.primID());
                 const Vec3fa v0 = mesh->vertex(tri.v[0]);
                 const Vec3fa v1 = mesh->vertex(tri.v[1]);
@@ -382,21 +382,21 @@ namespace embree
                progress,
 	       prims,pinfo,BVH8::N,BVH8::maxBuildDepthLeaf,sahBlockSize,minLeafSize,maxLeafSize,BVH8::travCost,intCost);
 	    bvh->set(root,pinfo.geomBounds,pinfo.size());
-            
+
 #if ROTATE_TREE
-            for (int i=0; i<ROTATE_TREE; i++) 
+            for (int i=0; i<ROTATE_TREE; i++)
               BVH8Rotate::rotate(bvh,bvh->root);
             bvh->clearBarrier(bvh->root);
 #endif
-            
+
             bvh->layoutLargeNodes(pinfo.size()*0.005f);
 
             if ((bvh->device->benchmark || bvh->device->verbosity(1)) && mesh == nullptr) dt = getSeconds()-t0;
-            
+
 #if PROFILE
             dt = timer.avg();
-        }); 
-#endif	
+        });
+#endif
 
 	/* clear temporary data for static geometry */
 	bool staticGeom = mesh ? mesh->isStatic() : scene->isStatic();