:segmentation: is ndarray with 0 and 1
    :voxelsize_mm: is array of three numbers specifiing size of voxel for each
        axis
    :aoi: is specify area of interest. It is ndarray with 0 and 1
    :sond_raster_mm: unimplemented. It is parametr of sonds design

def surface_density(segmentation, voxelsize_mm, aoi=None, sond_raster_mm=None):
    """
    :segmentation: is ndarray with 0 and 1
    :voxelsize_mm: is array of three numbers specifiing size of voxel for each
        axis
    :aoi: is specify area of interest. It is ndarray with 0 and 1
    :sond_raster_mm: unimplemented. It is parametr of sonds design
    """

    axis = 0
    if sond_raster_mm is None:
        sond_raster_mm = voxelsize_mm
    if aoi is None:
        aoi = np.ones(segmentation.shape)

    im_edg = find_edge(segmentation, axis=axis)
    im_edg = im_edg * aoi
    im_sond, aoi_sond = bufford_needle_sond(
        im_edg, voxelsize_mm, sond_raster_mm, axis=axis, aoi=aoi)

# isotropic fakir  - kubinova janacek
    # est S = 2 \frac{1}{n} \sum_{i=1}^{n} \frac{v}{l_i} \cdot l_i

# celkova delka sond
    # n_needle = (im_sond.shape[1] * im_sond.shape[2])
    # one_needle_l = im_sond.shape[0] * voxelsize_mm[0]
    # length = n_needle * one_needle_l
    length = np.sum(aoi_sond > 0) * voxelsize_mm[0]

# inverse of the probe per unit volume v/l_i
    # ippuv = (
    #     (np.prod(sond_raster_mm) * im_sond.shape[axis])
    #     /
    #     (sond_raster_mm[axis] * im_sond.shape[axis])
    # )
# Pocet pruseciku
    # Ii = np.sum(np.abs(im_sond))
    Ii = np.sum(np.abs(im_sond))

    # import sed3
    # ed = sed3.sed3(im_sond)
    # ed.show()

    # Kubinova2001
    # print "Ii = ", Ii
    Sv = 2.0 * Ii / length
    # import ipdb; ipdb.set_trace() #  noqa BREAKPOINT

    return Sv