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)