         // cout << "\n ** Before Leaf removal **\n" << ToString(true) << endl;
 
         // *** Remove leafs ***
-        CreateLeafNodes();        
-
+        CreateLeafNodes(rayMortonCodes);
+        
         // cout << "\n ** After Leaf removal **\n" << ToString(true) << endl;
     }
     CHECK_FOR_CUDA_ERROR();

    // Compute bounds of geometry partitions and do coarse level ray elimination
    // before intersection.

    // cout << "\n *** GORELESS VICTORY ***\n" << endl;
}


struct MortonIsNodeLeaf {
    // TODO take morton codes into account. If rayPartition.x and
    // rayPartition.y-1 points to the same morton code, then we can't possibly
    // subdivide further.

    __host__ __device__
    bool operator()(const uint2 rayPartition, const uint2 spherePartition) const {
        const float rayCount = (float)(rayPartition.y - rayPartition.x);
        const float sphereCount = (float)(spherePartition.y - spherePartition.x);
        
        //return rayCount * sphereCount <= 16.0f * (rayCount + sphereCount);
        return rayCount * sphereCount <= 11.0f * (rayCount + sphereCount);
    }
};


    out << i << ": [rays: [" << rayPartition.x << " -> " << rayPartition.y << 
        "], spheres: [" << spherePartition.x << " -> " << spherePartition.y << "]]";
    if (verbose) {
        Ray rMin = rays->GetLeafRay(rayPartition.x);
        Ray rMax = rays->GetLeafRay(rayPartition.y-1);
        
        RayMortonCoder rayMortonCoder(spheresGeom->GetBounds());
        MortonCode cMin = rayMortonCoder(rMin.origin, rMin.direction);
        MortonCode cMax = rayMortonCoder(rMin.origin, rMin.direction);
        out << "\n  Rays: [min: " << cMin << "\n         max: " << cMax << "]";
        out << "\n  Ray bound: " << MortonBound::LowestCommonBound(cMin, cMax);
        
        // out << "\n  Rays: [";
        // for (int r = rayPartition.x; r < rayPartition.y; ++r) {
        //     out << rays->GetLeafRay(r).id;
