docs: 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