                     ? min(d, get<axis>(next.max) - get<axis>(cell.max))
                     : max(d, get<axis>(next.min) - get<axis>(cell.min));
 
-                int first_ref = cell.begin;
-                for (int p = next.begin; p < next.end; p++) {
-                    auto ref = refs[p];
-                    auto found = bisection(refs + first_ref, cell.end - first_ref, ref);
-                    first_ref = found + 1 + first_ref;
-                    // If the reference is not in the cell we try to expand
-                    if (found < 0) {
-                        auto prim = load_prim(prims + ref);
-                        auto cur = dir ? max_bb : min_bb;
-                        auto left = d, right = dir ? 1 : -1;
-                        // Using bisection, find the offset by which we can overlap the neighbour
-                        while (dir ? (left >= right) : (left <= right)) {
-                            auto m = (left + right) / 2;
-                            if (axis == 0) cur.x = grid_min.x + cell_size.x * ((dir ? cell.max.x : cell.min.x) + m);
-                            if (axis == 1) cur.y = grid_min.y + cell_size.y * ((dir ? cell.max.y : cell.min.y) + m);
-                            if (axis == 2) cur.z = grid_min.z + cell_size.z * ((dir ? cell.max.z : cell.min.z) + m);
-                            if (intersect_prim_box(prim, dir ? BBox(min_bb, cur) : BBox(cur, max_bb))) {
-                                left = m + (dir ? -1 : 1);
-                            } else {
-                                right = m + (dir ? 1 : -1);
-                            }
-                        }
-                        d = left;
+                auto cell_bbox = BBox(grid_min + cell_size * vec3(cell.min),
+                                      grid_min + cell_size * vec3(cell.max));
+
+                int i = cell.begin, j = next.begin;
+                while (i < cell.end && j < next.end) {
+                    auto a = refs[i];
+                    auto b = refs[j];
+                    if (a < b) {
+                        // Present only in the current cell
+                        i++;
+                    } else if (a > b) {
+                        // Present only in the next cell
+                        d = compute_overlap<axis, dir>(load_prim(prims + b), cell, cell_bbox, d);
                         if (d == 0) break;
+                        j++;
+                    } else {
+                        // Present in both cells
+                        i++, j++;
                     }
                 }
+
+                while (j < next.end && d != 0)
+                    d = compute_overlap<axis, dir>(load_prim(prims + refs[j++]), cell, cell_bbox, d);
+
+                if (d == 0) break;
             }
-
-            if (d == 0) break;
 
             k1 = get<axis1>(next.max) - i;
             k2 = min(k2, get<axis2>(next.max) - j);

                        : max(d, get<axis>(next.min) - get<axis>(cell.min));
                } else {
                    d = 0;
                }
            } else {
                auto min_bb = grid_min + vec3(cell.min) * cell_size;

    return d;
}

template <int axis, bool last_iter, typename Primitive>
__global__ void overlap_step(const Entry* __restrict__ entries,
                             const int* __restrict__ refs,
                             const Primitive* prims,
