docs: def __surface_distances(result, reference, voxelspacing=None, connectivity=1):
      """
      The distances between the surface voxel of binary objects in result and their
      nearest partner surface voxel of a binary object in reference.
      """
      result = numpy.atleast_1d(result.astype(numpy.bool))
      reference = numpy.atleast_1d(reference.astype(numpy.bool))
      if voxelspacing is not None:
          voxelspacing = _ni_support._normalize_sequence(voxelspacing, result.ndim)
          voxelspacing = numpy.asarray(voxelspacing, dtype=numpy.float64)
          if not voxelspacing.flags.contiguous:
              voxelspacing = voxelspacing.copy()
              
      # binary structure
      footprint = generate_binary_structure(result.ndim, connectivity)
      
      # test for emptiness
      if 0 == numpy.count_nonzero(result): 
          raise RuntimeError('The first supplied array does not contain any binary object.')
      if 0 == numpy.count_nonzero(reference): 
          raise RuntimeError('The second supplied array does not contain any binary object.')    
              
      # extract only 1-pixel border line of objects
      result_border = result ^ binary_erosion(result, structure=footprint, iterations=1)
      reference_border = reference ^ binary_erosion(reference, structure=footprint, iterations=1)
      
      # compute average surface distance        
      # Note: scipys distance transform is calculated only inside the borders of the
      #       foreground objects, therefore the input has to be reversed
      dt = distance_transform_edt(~reference_border, sampling=voxelspacing)
      sds = dt[result_border]
      
      return sds