docs: def boundary_stawiaski(graph, label_image, gradient_image): # label image is not required to hold continuous ids or to start from 1
      r"""
      Boundary term based on the sum of border voxel pairs differences.
       
      An implementation of the boundary term in [1]_, suitable to be used with the `~medpy.graphcut.generate.graph_from_labels` function.
      
      Determines for each two supplied regions the voxels forming their border assuming
      :math:`ndim*2`-connectedness (e.g. :math:`3*2=6` for 3D). From the gradient magnitude values of each
      end-point voxel the border-voxel pairs, the highest one is selected and passed to a
      strictly positive and decreasing function :math:`g(x)`, which is defined as:
      
      .. math::
      
          g(x) = \left(\frac{1}{1+|x|}\right)^k
      
      ,where :math:`k=2`. The final weight :math:`w_{i,j}` between two regions :math:`r_i` and
      :math:`r_j` is then determined by the sum of all these neighbour values:
      
      .. math::
      
          w_{i,j} = \sum_{e_{m,n}\in F_{(r_i,r_j)}}g(\max(|I(m)|,|I(n)|))
      
      , where :math:`F_{(r_i,r_j)}` is the set of border voxel-pairs :math:`e_{m,n}` between
      the regions :math:`r_i` and :math:`r_j` and :math:`|I(p)|` the absolute of the gradient
      magnitude at the voxel :math:`p`
      
      This boundary_function works as an edge indicator in the original image. In simpler
      words the weight (and therefore the energy) is obtained by summing the local contrast
      along the boundaries between two regions.
      
      Parameters
      ----------
      graph : GCGraph
          The graph to add the weights to.
      label_image : ndarray
          The label image. Must contain consecutively labelled regions starting from index 1.
      gradient_image : ndarray
          The gradient image.
      
      Notes
      -----
      This function requires the gradient magnitude image of the original image to be passed
      along. That means that `~medpy.graphcut.generate.graph_from_labels` has to be called
      with ``boundary_term_args`` set to the gradient image. This can be obtained e.g. with
      `generic_gradient_magnitude` and `prewitt` from `scipy.ndimage`.
      
      This function is tested on 2D and 3D images and theoretically works for all dimensionalities. 
      
      References
      ----------
      .. [1] Stawiaski J., Decenciere E., Bidlaut F. "Interactive Liver Tumor Segmentation
             Using Graph-cuts and watershed" MICCAI 2008 participation
      """
      # convert to arrays if necessary
      label_image = scipy.asarray(label_image)
      gradient_image = scipy.asarray(gradient_image)
      
      if label_image.flags['F_CONTIGUOUS']: # strangely, this one is required to be ctype ordering
          label_image = scipy.ascontiguousarray(label_image)
          
      __check_label_image(label_image)
          
      for dim in range(label_image.ndim):
          # prepare slicer for all minus last and all minus first "row"
          slicer_from = [slice(None)] * label_image.ndim
          slicer_to = [slice(None)] * label_image.ndim
          slicer_from[dim] = slice(None, -1)
          slicer_to[dim] = slice(1, None)
          # slice views of keys
          keys_from = label_image[slicer_from]
          keys_to = label_image[slicer_to]
          # determine not equal keys
          valid_edges = keys_from != keys_to
          # determine largest gradient
          gradient_max = numpy.maximum(numpy.abs(gradient_image[slicer_from]), numpy.abs(gradient_image[slicer_to]))[valid_edges]
          # determine key order
          keys_max = numpy.maximum(keys_from, keys_to)[valid_edges]
          keys_min = numpy.minimum(keys_from, keys_to)[valid_edges]
          # set edges / nweights
          for k1, k2, val in zip(keys_min, keys_max, gradient_max):
              weight = math.pow(1./(1. + val), 2) # weight contribution of a single pixel
              weight = max(weight, sys.float_info.min)
              graph.set_nweight(k1 - 1 , k2 - 1, weight, weight)