     }
 
     __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) {
+                                       const size_t boundsSize,
+                                       MortonCode* sharedMin) {
     
         __shared__ volatile unsigned int mins[128];
         __shared__ volatile unsigned int maxs[128];
 
-        if (threadIdx.x < intermediateBoundsSize) {
-            MortonBound b = intermediateBounds[threadIdx.x];
-            mins[threadIdx.x] = b.min;
-            maxs[threadIdx.x] = b.max;
-        } else
-            mins[threadIdx.x] = maxs[threadIdx.x] = MortonCode::EncodeAxis(6);
-        __syncthreads();
+        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];

        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);

    }


        }
    }


    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)
