     leafNodeIndices[0] = 0;
     thrust::transform_inclusive_scan(isLeaf.begin(), isLeaf.end(), leafNodeIndices.begin()+1, 
                                      BoolToInt(), thrust::plus<unsigned int>());
-    cout << "leafNodeIndices:\n" << leafNodeIndices << endl;
+    cout << "leafNodeEntries:\n" << leafNodeIndices << endl;
     const size_t nNewLeafs = leafNodeIndices[leafNodeIndices.size()-1];
     cout << "nNewLeafs: " << nNewLeafs << endl;
     
 
     if (nNewLeafs == 0) { 
-        cout << "No new leafs" << endl;
+        // cout << "No new leafs" << endl;
         return false;
     }
 

        bounds.resize(nextActiveNodes);
        thrust::transform(nextRayPartitions.begin() + leafNodes, nextRayPartitions.end(), 
                          bounds.begin(), CreateBoundsFromPartitions(rayMortonCodes));
        std::cout << "bounds:\n" << bounds << std::endl;

        static HyperCubes hyperCubes(128); hyperCubes.Resize(bounds.size());
        thrust::transform(bounds.begin(), bounds.end(), 

    spherePartitions.swap(nextSpherePartitions);
    rayPartitions.swap(nextRayPartitions);

    // cout << "new spherePartitions:\n" << spherePartitions << endl;
    // cout << "new rayPartitions:\n" << rayPartitions << endl;
    

    leafSphereIndices += nNewLeafIndices;
    leafNodes += nNewLeafs;
    // cout << "new leaf count: " << leafNodes << endl;
    // cout << "leafSphereIndices:\n" << leafSphereIndices << endl;

    return true;
}


__global__
void ExhaustiveIntersection(const float4* const rayOrigins, 
                            float4* rayDirs,
                            const unsigned int* const nodes,
                            const uint2* const spherePartitionings,
                            const unsigned int* const sphereIndices,
                            const Sphere* const spheres,
                            const unsigned int nRays,
                            unsigned int *hitIDs) {
    
    unsigned int id = threadIdx.x + blockDim.x * blockIdx.x;
    if (id >= nRays) return;

    const unsigned int node = nodes[id];

    const float3 origin = make_float3(rayOrigins[id]);
    const float3 dir = make_float3(rayDirs[id]);
        
    const uint2 spherePartitioning = spherePartitionings[node];
    float hitT = 1e30f;
    unsigned int hitID = SpheresGeometry::MISSED;
    
    for (unsigned int g = spherePartitioning.x; g < spherePartitioning.y; ++g) {
        const unsigned int sphereId = sphereIndices[g];
        const Sphere s = spheres[sphereId];
        const float t = s.Intersect(origin, dir);
        if (0 < t && t < hitT) {
            hitID = sphereId;
            hitT = t;
        }
    }
    
    rayDirs[id] = make_float4(dir, hitT);
    hitIDs[id] = hitID;
}

void MortonDacrtNodes::FindIntersections(thrust::device_vector<unsigned int>& hitIDs) {
    const unsigned int nActiveRays = rays->LeafRays();
    hitIDs.resize(nActiveRays);
    
    // Compute which nodes each ray belongs to. This should be done using work
    // queues on high end cards.
    static thrust::device_vector<unsigned int> owners(nActiveRays);
    owners.resize(nActiveRays);
    DacrtNodes::CalcOwners(RayPartitionsBegin(), RayPartitionsEnd(), owners);

    // Perform exhaustive intersection.
    struct cudaFuncAttributes funcAttr;
    cudaFuncGetAttributes(&funcAttr, ExhaustiveIntersection);
    unsigned int blocksize = funcAttr.maxThreadsPerBlock > 256 ? 256 : funcAttr.maxThreadsPerBlock;
    unsigned int blocks = (nActiveRays / blocksize) + 1;
    ExhaustiveIntersection<<<blocks, blocksize>>>
        (RawPointer(Rays::GetOrigins(rays->BeginLeafRays())),
         RawPointer(Rays::GetDirections(rays->BeginLeafRays())),
         RawPointer(owners),
         RawPointer(spherePartitions),
         RawPointer(sphereIndices),
         RawPointer(spheresGeom->spheres),
         nActiveRays,
         RawPointer(hitIDs));
    CHECK_FOR_CUDA_ERROR();
}


