code: def shape_based_slice_insertation(sl1, sl2, dim, nslices, order=3):
      """
      Insert `nslices` new slices between `sl1` and `sl2` along dimension `dim` using shape
      based binary interpolation.
  
      Extrapolation is handled adding `nslices`/2 step-wise eroded copies of the last slice
      in each direction.
  
      Parameters
      ----------
      sl1 : array_like
          First slice. Treated as binary data.
      sl2 : array_like
          Second slice. Treated as binary data.
      dim : int
          The new dimension along which to add the new slices.
      nslices : int
          The number of slices to add.
      order : int
          The b-spline interpolation order for re-sampling the distance maps.
  
      Returns
      -------
      out : ndarray
          A binary image of size `sl1`.shape() extend by `nslices`+2 along the new
          dimension `dim`. The border slices are the original slices `sl1` and `sl2`.
      """
      sl1 = sl1.astype(numpy.bool_)
      sl2 = sl2.astype(numpy.bool_)
  
      # extrapolation through erosion
      if 0 == numpy.count_nonzero(sl1):
          slices = [sl1]
          for _ in range(nslices / 2):
              slices.append(numpy.zeros_like(sl1))
          for i in range(1, nslices / 2 + nslices % 2 + 1)[::-1]:
              slices.append(binary_erosion(sl2, iterations=i))
          slices.append(sl2)
          return numpy.rollaxis(numpy.asarray(slices), 0, dim + 1)
          # return numpy.asarray([sl.T for sl in slices]).T
      elif 0 == numpy.count_nonzero(sl2):
          slices = [sl1]
          for i in range(1, nslices / 2 + 1):
              slices.append(binary_erosion(sl1, iterations=i))
          for _ in range(0, nslices / 2 + nslices % 2):
              slices.append(numpy.zeros_like(sl2))
          slices.append(sl2)
          return numpy.rollaxis(numpy.asarray(slices), 0, dim + 1)
          # return numpy.asarray([sl.T for sl in slices]).T
  
      # interpolation shape based
      # note: distance_transform_edt shows strange behaviour for ones-arrays
      dt1 = distance_transform_edt(~sl1) - distance_transform_edt(sl1)
      dt2 = distance_transform_edt(~sl2) - distance_transform_edt(sl2)
  
      slicer = [slice(None)] * dt1.ndim
      slicer = slicer[:dim] + [numpy.newaxis] + slicer[dim:]
      out = numpy.concatenate((dt1[slicer], dt2[slicer]), axis=dim)
      zoom_factors = [1] * dt1.ndim
      zoom_factors = (
          zoom_factors[:dim]
          + [(nslices + 2) / 2.0]
          + zoom_factors[dim:]
      )
      out = zoom(out, zoom_factors, order=order)
  
      return out <= 0
  
file_path: medpy/bin/medpy_binary_resampling.py
imports: import os
  import logging
  import argparse
  import numpy
  from scipy.ndimage import zoom
  from scipy.ndimage import distance_transform_edt, binary_erosion
  from scipy.ndimage import label
  from medpy.core import Logger
  from medpy.filter import resample, bounding_box
  from medpy.utilities import argparseu
  from medpy.io import load, save, header
label_desc: Use image node_id 4 for calling a global function with example usage: shape_based_slice_insertation(sl1, sl2, dim, nslices, order) and returns: unknown, numpy, numpy