docs: 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])