docs: def boundary_difference_of_means(graph, label_image, original_image): # label image is not required to hold continuous ids or to start from 1
      r"""
      Boundary term based on the difference of means between adjacent image regions.
      
      An implementation of the boundary term, suitable to be used with the `~medpy.graphcut.generate.graph_from_labels` function.
      
      This simple energy function computes the mean values for all regions. The weights of
      the edges are then determined by the difference in mean values.
      
      The graph weights generated have to be strictly positive and preferably in the
      interval :math:`(0, 1]`. To ensure this, the maximum possible difference in mean values is computed as:
      
      .. math::
      
          \alpha = \|\max \bar{I} - \min \bar{I}\|
      
      , where :math:`\min \bar{I}` constitutes the lowest mean intensity value of all regions in
      the image, while :math:`\max \bar{I}` constitutes the highest mean intensity value With this
      value the weights between a region :math:`x` and its neighbour :math:`y` can be computed:
      
      .. math::
      
          w(x,y) = \max \left( 1 - \frac{\|\bar{I}_x - \bar{I}_y\|}{\alpha}, \epsilon \right)
      
      where :math:`\epsilon` is the smallest floating point step and thus :math:`w(x,y) \in (0, 1]` holds true.
      
      Parameters
      ----------
      graph : GCGraph
          The graph to add the weights to.
      label_image : ndarray
          The label image.
      original_image : ndarray
          The original image.
      
      Notes
      -----
      This function requires 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
      original image. 
      
      This function is tested on 2D and 3D images and theoretically works for all dimensionalities. 
      """
      # convert to arrays if necessary
      label_image = scipy.asarray(label_image)
      original_image = scipy.asarray(original_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)
      
      # create a lookup-table that translates from a label id to its position in the sorted unique vector
      labels_unique = scipy.unique(label_image)
      
      # compute the mean intensities of all regions
      # Note: Bug in mean implementation: means over labels is only computed if the indexes are also supplied
      means = scipy.ndimage.measurements.mean(original_image, labels=label_image, index=labels_unique)
      
      # compute the maximum possible intensity difference
      max_difference = float(abs(min(means) - max(means)))
  
      # create a lookup table that relates region ids to their respective intensity values
      means = dict(list(zip(labels_unique, means)))
  
      # get the adjuncancy of the labels
      edges = __compute_edges(label_image)
      
      # compute the difference of means for each adjunct region and add it as a tuple to the dictionary
      if 0. == max_difference: # special case when the divider is zero and therefore all values can be assured to equal zero
          for edge in edges:
              graph.set_nweight(edge[0] - 1, edge[1] - 1, sys.float_info.min, sys.float_info.min)
      else:    
          # compute the difference of means for each adjunct region and add it as a tuple to the dictionary
          for edge in edges:
              value = max(1. - abs(means[edge[0]] - means[edge[1]]) / max_difference, sys.float_info.min)
              graph.set_nweight(edge[0] - 1, edge[1] - 1, value, value)