def boundary_stawiaski_directed(graph, label_image, xxx_todo_changeme): # 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, directed version.
    
    An implementation of the boundary term in [1]_, suitable to be used with the
    `~medpy.graphcut.generate.graph_from_labels` function.
    
    The basic definition of this term is the same as for `boundary_stawiaski`, but the
    edges of the created graph will be directed.
    
    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.
    
    When the ``directedness`` parameter is set to zero, the resulting graph will be undirected
    and the behaviour equals `boundary_stawiaski`.
    When it is set to a positive value, light-to-dark transitions are favored i.e. voxels
    with a lower intensity (darker) than the objects tend to be assigned to the object.
    The boundary term is thus changed to:
    
    .. math::
    
          g_{ltd}(x) = \left\{
              \begin{array}{l l}
                g(x) + \beta & \quad \textrm{if $I_i > I_j$}\\
                g(x) & \quad \textrm{if $I_i \leq I_j$}\\
              \end{array} \right.

    With a negative value for ``directedness``, the opposite effect can be achieved i.e.
    voxels with a higher intensity (lighter) than the objects tend to be assigned to the
    object. The boundary term is thus changed to
    
    .. math::
    
      g_{dtl} = \left\{
          \begin{array}{l l}
            g(x) & \quad \textrm{if $I_i > I_j$}\\
            g(x) + \beta & \quad \textrm{if $I_i \leq I_j$}\\
          \end{array} \right.

    Subsequently the :math:`g(x)` in the computation of :math:`w_{i,j}` is substituted by
    :math:`g_{ltd}` resp. :math:`g_{dtl}`. The value :math:`\beta` determines the power of the
    directedness and corresponds to the absolute value of the supplied ``directedness``
    parameter. Experiments showed values between 0.0001 and 0.0003 to be good candidates.
    
    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.
    directedness : integer
        The weight of the directedness, a positive number to favour
        light-to-dark and a negative to dark-to-light transitions. See function
        description for more details.
    
    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    
    """
    (gradient_image, directedness) = xxx_todo_changeme
    label_image = scipy.asarray(label_image)
    gradient_image = scipy.asarray(gradient_image)
    
    if label_image.flags['F_CONTIGUOUS']: # strangely one this one is required to be ctype ordering
        label_image = scipy.ascontiguousarray(label_image)
        
    __check_label_image(label_image)
        
    beta = abs(directedness)
        
    def addition_directed_ltd(key1, key2, v1, v2, dic): # for light-to-dark # tested
        "Takes a key defined by two uints, two voxel intensities and a dict to which it adds g(v1, v2)."
        if not key1 == key2: # do not process voxel pairs which belong to the same region
            # The function used to compute the weight contribution of each voxel pair
            weight = math.pow(1./(1. + max(abs(v1), abs(v2))), 2)
            # ensure that no value is zero; this can occur due to rounding errors
            weight = max(weight, sys.float_info.min)
            # add weighted values to already existing edge
            if v1 > v2: graph.set_nweight(key1 - 1, key2 - 1, min(1, weight + beta), weight)
            else: graph.set_nweight(key1 - 1, key2 - 1, weight, min(1, weight + beta))
            
    def addition_directed_dtl(key1, key2, v1, v2): # for dark-to-light # tested
        "Takes a key defined by two uints, two voxel intensities and a dict to which it adds g(v1, v2)."
        if not key1 == key2: # do not process voxel pairs which belong to the same region
            # The function used to compute the weight contribution of each voxel pair
            weight = math.pow(1./(1. + max(abs(v1), abs(v2))), 2)
            # ensure that no value is zero; this can occur due to rounding errors
            weight = max(weight, sys.float_info.min)
            # add weighted values to already existing edge
            if v1 > v2: graph.set_nweight(key1 - 1, key2 - 1, weight, min(1, weight + beta))
            else: graph.set_nweight(key1 - 1, key2 - 1, min(1, weight + beta), weight)
                                                  
    # pick and vectorize the function to achieve a speedup
    if 0 > directedness:
        vaddition = scipy.vectorize(addition_directed_dtl)
    else:
        vaddition = scipy.vectorize(addition_directed_ltd)
    
    # iterate over each dimension
    for dim in range(label_image.ndim):
        slices_x = []
        slices_y = []
        for di in range(label_image.ndim):
            slices_x.append(slice(None, -1 if di == dim else None))
            slices_y.append(slice(1 if di == dim else None, None))
        vaddition(label_image[slices_x],
                  label_image[slices_y],
                  gradient_image[slices_x],
                  gradient_image[slices_y])
r"""
    Boundary term based on the sum of border voxel pairs differences, directed version.
    
    An implementation of the boundary term in [1]_, suitable to be used with the
    `~medpy.graphcut.generate.graph_from_labels` function.
    
    The basic definition of this term is the same as for `boundary_stawiaski`, but the
    edges of the created graph will be directed.
    
    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.
    
    When the ``directedness`` parameter is set to zero, the resulting graph will be undirected
    and the behaviour equals `boundary_stawiaski`.
    When it is set to a positive value, light-to-dark transitions are favored i.e. voxels
    with a lower intensity (darker) than the objects tend to be assigned to the object.
    The boundary term is thus changed to:
    
    .. math::
    
          g_{ltd}(x) = \left\{
              \begin{array}{l l}
                g(x) + \beta & \quad \textrm{if $I_i > I_j$}\\
                g(x) & \quad \textrm{if $I_i \leq I_j$}\\
              \end{array} \right.

    With a negative value for ``directedness``, the opposite effect can be achieved i.e.
    voxels with a higher intensity (lighter) than the objects tend to be assigned to the
    object. The boundary term is thus changed to
    
    .. math::
    
      g_{dtl} = \left\{
          \begin{array}{l l}
            g(x) & \quad \textrm{if $I_i > I_j$}\\
            g(x) + \beta & \quad \textrm{if $I_i \leq I_j$}\\
          \end{array} \right.

    Subsequently the :math:`g(x)` in the computation of :math:`w_{i,j}` is substituted by
    :math:`g_{ltd}` resp. :math:`g_{dtl}`. The value :math:`\beta` determines the power of the
    directedness and corresponds to the absolute value of the supplied ``directedness``
    parameter. Experiments showed values between 0.0001 and 0.0003 to be good candidates.
    
    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.
    directedness : integer
        The weight of the directedness, a positive number to favour
        light-to-dark and a negative to dark-to-light transitions. See function
        description for more details.
    
    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