def depth_to_voxels(
    ray_voxel_count,
    ray_voxel_indices,
    rays_idxs,
    voxel_grid,
    points,
    S,
    S_new,
    single_ray_depth_to_voxels,
    gamma=None
):
    """Compute the depth probability of each voxel based on S the probability
    distribution on points.

    Arguments
    ---------
        ray_voxel_count: array(shape=(N,), int), The number of voxels
                         intersected by each ray
        ray_voxel_indices: array(shape=(N, M, 3), int), The indices in the
                           voxel grid per ray
        rays_idxs: array(shape=(N,), int), The indices of the valid rays
        voxel_gird: array(shape=(3, D1, D2, D3)), The coordinates of the
                    centers of all voxels in the voxel grid of size
                    (D1, D2, D3)
        points: array(shape=(4, N, D), float32), D points sampled on each of
                the N rays
        S: array(shape=(N, D), float32), A depth probability distribution for
           each of the N rays
        single_ray_depth_to_voxels: callable that sets the depth to voxels
                                    sampling method
        gamma: int, gamma prior for the gaussian distribution, when we do the
               mapping using the kernel density estimation
    """
    # Extract the numbers N, M, D
    N, M, _ = ray_voxel_indices.shape
    _, _, D = points.shape

    # Compute the voxel probabilities
    S_new.fill(0)
    for r, ray_idx in enumerate(rays_idxs):
        C = ray_voxel_count[r]
        idxs = ray_voxel_indices[r, :C]
        S_new[r, :C] = single_ray_depth_to_voxels(
            voxel_grid[:, idxs[:, 0], idxs[:, 1], idxs[:, 2]],
            points[:3, ray_idx],
            S[ray_idx],
            gamma=gamma
        )

    return S_new