<SYSTEM_TASK:>
Solve the following problem using Python, implementing the functions described below, one line at a time
<END_TASK>
<USER_TASK:>
Description:
 def overlapping_spheres(shape: List[int], radius: int, porosity: float, iter_max: int = 10, tol: float = 0.01):
 r""" Generate a packing of overlapping mono-disperse spheres Parameters shape : list The size of the image to generate in [Nx, Ny, Nz] where Ni is the number of voxels in the i-th direction. radius : scalar The radius of spheres in the packing. porosity : scalar The porosity of the final image, accurate to the given tolerance. iter_max : int Maximum number of iterations for the iterative algorithm that improves the porosity of the final image to match the given value. tol : float Tolerance for porosity of the final image compared to the given value. Returns ------- image : ND-array A boolean array with ``True`` values denoting the pore space Notes ----- This method can also be used to generate a dispersion of hollows by treating ``porosity`` as solid volume fraction and inverting the returned image. """

    shape = sp.array(shape)
    if sp.size(shape) == 1:
        shape = sp.full((3, ), int(shape))
    ndim = (shape != 1).sum()
    s_vol = ps_disk(radius).sum() if ndim == 2 else ps_ball(radius).sum()

    bulk_vol = sp.prod(shape)
    N = int(sp.ceil((1 - porosity)*bulk_vol/s_vol))
    im = sp.random.random(size=shape)

    # Helper functions for calculating porosity: phi = g(f(N))
    f = lambda N: spim.distance_transform_edt(im > N/bulk_vol) < radius
    g = lambda im: 1 - im.sum() / sp.prod(shape)

    # # Newton's method for getting image porosity match the given
    # w = 1.0                         # Damping factor
    # dN = 5 if ndim == 2 else 25     # Perturbation
    # for i in range(iter_max):
    #     err = g(f(N)) - porosity
    #     d_err = (g(f(N+dN)) - g(f(N))) / dN
    #     if d_err == 0:
    #         break
    #     if abs(err) <= tol:
    #         break
    #     N2 = N - int(err/d_err)   # xnew = xold - f/df
    #     N = w * N2 + (1-w) * N

    # Bisection search: N is always undershoot (bc. of overlaps)
    N_low, N_high = N, 4*N
    for i in range(iter_max):
        N = sp.mean([N_high, N_low], dtype=int)
        err = g(f(N)) - porosity
        if err > 0:
            N_low = N
        else:
            N_high = N
        if abs(err) <= tol:
            break

    return ~f(N)