def _ul_subvoxel_overlap(xs, x1, x2):
    """For an interval [x1, x2], return the index of the upper limit of the
    overlapping subvoxels whose borders are defined by the elements of xs."""
    xmax = max(x1, x2)
    if xmax >= xs[-1]:
        return len(xs) - 1
    elif xmax <= xs[0]:
        ul = 0
        return ul
    else:
        for i, x in enumerate(xs):
            if not x < xmax:
                ul = i
                return ul
        ul = len(xs) - 1
        return ul