         }
     }
 
-    __global__
-    void ReduceMinMaxMortonByAxisPass1(const unsigned int* const mortonCodes,
-                                       const size_t inputSize,
-                                       MortonBound* intermediateBounds,
-                                       MortonBound* bounds,
-                                       const size_t boundsSize) {
-
-        // Fill bounds with default values
-        for (unsigned int i = blockDim.x * blockIdx.x + threadIdx.x; i < boundsSize; i += blockDim.x * gridDim.x)
-            bounds[i] = MortonBound::Create(0xDFFFFFFF, 0xC0000000);
-
-        __shared__ volatile unsigned int min[256];
-        __shared__ volatile unsigned int max[256];
-    
-        const size_t beginIndex = inputSize * blockIdx.x / gridDim.x;
-        const size_t endIndex = inputSize * (blockIdx.x + 1) / gridDim.x;
-    
-        size_t currentIndex = beginIndex + threadIdx.x;
-
-        // Fill initial values
-        min[threadIdx.x] = max[threadIdx.x] = currentIndex < endIndex ? MortonCode(mortonCodes[currentIndex]) : MortonCode::EncodeAxis(6);
-        currentIndex += blockDim.x;
-        __syncthreads();
-    
-        // While still values left, load them and perform reduction
-        while (currentIndex < endIndex) {
-            // Fetch new data from global mem to shared mem
-            unsigned int code = currentIndex < endIndex ? MortonCode(mortonCodes[currentIndex]) : MortonCode::EncodeAxis(6);
-            min[blockDim.x + threadIdx.x] = max[blockDim.x + threadIdx.x] = code;
-            __syncthreads();
-
-            PerformMinMaxMortonByAxisReduction(threadIdx.x, min, max, bounds);
-        }
-        __syncthreads();
-
-        // Reduce 128 values left in shared memory
-        if (threadIdx.x >= 64) return;
-        PerformMinMaxMortonByAxisReduction(threadIdx.x, min, max, bounds);
-
-        if (threadIdx.x >= 32) return;
-        PerformMinMaxMortonByAxisReduction(threadIdx.x, min, max, bounds);
-
-        if (threadIdx.x >= 16) return;
-        PerformMinMaxMortonByAxisReduction(threadIdx.x, min, max, bounds);
-
-        if (threadIdx.x >= 8) return;
-        PerformMinMaxMortonByAxisReduction(threadIdx.x, min, max, bounds);
-
-        if (threadIdx.x >= 4) return;
-        PerformMinMaxMortonByAxisReduction(threadIdx.x, min, max, bounds);
-    
-        if (threadIdx.x >= 2) return;
-        // Store the last 2 values, since they may overlap with the first and last
-        // value in neighbouring blocks.
-        intermediateBounds[2 * blockIdx.x + threadIdx.x] = MortonBound::Create((unsigned int)min[threadIdx.x], 
-                                                                               (unsigned int)max[threadIdx.x]);
-    }
-
 
     void ReduceMinMaxMortonByAxis(thrust::device_vector<unsigned int>::iterator mortonBegin,
                                   thrust::device_vector<unsigned int>::iterator mortonEnd,
                                   thrust::device_vector<MortonBound>::iterator boundsBegin,
                                   thrust::device_vector<MortonBound>::iterator boundsEnd) {
+
+        /*
+        MortonCode code = (unsigned int)*mortonBegin;
+        MortonBound bound = MortonBound::Create(code, code);
+        ++mortonBegin;
+        while (mortonBegin != mortonEnd) {
+            code = *mortonBegin;
+            SignedAxis axis = bound.min.GetAxis();
+            if (axis == code.GetAxis()) {
+                bound.min = Morton::MinBy4(bound.min.WithoutAxis(), code.WithoutAxis()) + axis;
+                bound.max = Morton::MaxBy4(bound.max.WithoutAxis(), code.WithoutAxis()) + axis;
+            } else {
+                *(boundsBegin+(int)axis) = bound;
+                
+                bound = MortonBound::Create(code, code);
+            }
+            
+            ++mortonBegin;
+        }
+        SignedAxis axis = bound.min.GetAxis();
+        *(boundsBegin+(int)axis) = bound;
+
+        return;
+        */
 
         // Verify that CUDA is initialized
         if (!Meta::CUDA::initialized)

        const unsigned int blocks = Meta::CUDA::activeCudaDevice.multiProcessorCount;

        static thrust::device_vector<MortonBound> intermediateBounds(boundsSize * 2);
    
        ReduceMinMaxMortonByAxisPass1<<<blocks, blockDim>>>(RawPointer(mortonBegin), inputSize,
                                                            RawPointer(intermediateBounds),
                                                            thrust::raw_pointer_cast(&*boundsBegin), boundsSize);

        // std::cout << "intermediateBounds\n" << intermediateBounds << std::endl;
    
        // std::cout << "bounds:\n";
        // for (int i = 0; i < boundsSize; ++i) {
        //     std::cout << i << ": " << boundsBegin[i];
        //     if (i < boundsSize-1)
        //         std::cout << "\n";
        // }
        // std::cout << std::endl;
    
        ReduceMinMaxMortonByAxisPass2<<<1, 128>>>(RawPointer(intermediateBounds), intermediateBounds.size(),
                                                  thrust::raw_pointer_cast(&*boundsBegin), boundsSize);

    }


    }

    __global__
    void ReduceMinMaxMortonByAxisPass1(const unsigned int* const mortonCodes,
                                       const size_t inputSize,
                                       MortonBound* intermediateBounds,
                                       MortonBound* bounds,
                                       const size_t boundsSize) {

        // Fill bounds with default values
        for (unsigned int i = blockDim.x * blockIdx.x + threadIdx.x; i < boundsSize; i += blockDim.x * gridDim.x)
            bounds[i] = MortonBound::Create(0xDFFFFFFF, 0xC0000000);

        __shared__ volatile unsigned int min[256];
        __shared__ volatile unsigned int max[256];
    
        const size_t beginIndex = inputSize * blockIdx.x / gridDim.x;
        const size_t endIndex = inputSize * (blockIdx.x + 1) / gridDim.x;
    
        size_t currentIndex = beginIndex + threadIdx.x;

        // Fill initial values
        unsigned int lookupIndex = Min(currentIndex, endIndex-1); // Ensures that we pad with the last element
        min[threadIdx.x] = max[threadIdx.x] = MortonCode(mortonCodes[lookupIndex]);
        currentIndex += blockDim.x;
        __syncthreads();
    
        // While still values left, load them and perform reduction
        while (currentIndex < endIndex) {
            // Fetch new data from global mem to shared mem
            lookupIndex = Min(currentIndex, endIndex-1); // Ensures that we pad with the last element
            min[blockDim.x + threadIdx.x] = max[blockDim.x + threadIdx.x] = MortonCode(mortonCodes[lookupIndex]);
            __syncthreads();

            PerformMinMaxMortonByAxisReduction(threadIdx.x, min, max, bounds);

            currentIndex += blockDim.x;
        }
        __syncthreads();

        // Reduce 128 values left in shared memory
        if (threadIdx.x >= 64) return;
        PerformMinMaxMortonByAxisReduction(threadIdx.x, min, max, bounds);
        __syncthreads();

        if (threadIdx.x >= 32) return;
        PerformMinMaxMortonByAxisReduction(threadIdx.x, min, max, bounds);

        if (threadIdx.x >= 16) return;
        PerformMinMaxMortonByAxisReduction(threadIdx.x, min, max, bounds);

        if (threadIdx.x >= 8) return;
        PerformMinMaxMortonByAxisReduction(threadIdx.x, min, max, bounds);

        if (threadIdx.x >= 4) return;
        PerformMinMaxMortonByAxisReduction(threadIdx.x, min, max, bounds);
    
        if (threadIdx.x >= 2) return;
        // Store the last 2 values, since they may overlap with the first and last
        // value in neighbouring blocks.
        intermediateBounds[2 * blockIdx.x + threadIdx.x] = MortonBound::Create((unsigned int)min[threadIdx.x], 
                                                                               (unsigned int)max[threadIdx.x]);
    }

    __global__
    void ReduceMinMaxMortonByAxisPass2(MortonBound* intermediateBounds,
                                       const size_t intermediateBoundsSize,
                                       MortonBound* bounds,
                                       const size_t boundsSize,
                                       MortonCode* sharedMin) {
    
        __shared__ volatile unsigned int mins[128];
        __shared__ volatile unsigned int maxs[128];

        unsigned int index = Min((size_t)threadIdx.x, intermediateBoundsSize-1);
        MortonBound b = intermediateBounds[index];
        mins[threadIdx.x] = b.min;
        maxs[threadIdx.x] = b.max;
        __syncthreads();

        sharedMin[threadIdx.x] = mins[threadIdx.x];

        // Reduce 128 values left in shared memory
        if (threadIdx.x >= 64) return;
        PerformMinMaxMortonByAxisReduction(threadIdx.x, mins, maxs, bounds);
        __syncthreads();

        sharedMin[128 + threadIdx.x] = mins[threadIdx.x];
        sharedMin[128 + 64 + threadIdx.x] = mins[64 + threadIdx.x];

        if (threadIdx.x >= 32) return;
        PerformMinMaxMortonByAxisReduction(threadIdx.x, mins, maxs, bounds);

        for (unsigned int i = 0; i < 128; i+=32)
            sharedMin[256 + i + threadIdx.x] = mins[i + threadIdx.x];

        if (threadIdx.x >= 16) return;
        PerformMinMaxMortonByAxisReduction(threadIdx.x, mins, maxs, bounds);

        for (unsigned int i = 0; i < 128; i+=16)
            sharedMin[384 + i + threadIdx.x] = mins[i + threadIdx.x];

        if (threadIdx.x >= 8) return;
        PerformMinMaxMortonByAxisReduction(threadIdx.x, mins, maxs, bounds);
        for (unsigned int i = 0; i < 128; i+=8)
            sharedMin[512 + i + threadIdx.x] = mins[i + threadIdx.x];

        if (threadIdx.x >= 4) return;
        PerformMinMaxMortonByAxisReduction(threadIdx.x, mins, maxs, bounds);
        for (unsigned int i = 0; i < 128; i+=4)
            sharedMin[640 + i + threadIdx.x] = mins[i + threadIdx.x];
    
        if (threadIdx.x >= 2) return;
        PerformMinMaxMortonByAxisReduction(threadIdx.x, mins, maxs, bounds);
        for (unsigned int i = 0; i < 128; i+=2)
            sharedMin[768 + i + threadIdx.x] = mins[i + threadIdx.x];
    
        // Reduce the last value, comparing it to what is already stored in bounds
        const unsigned int lhsMin = mins[0];
