docs: def directionality(image, min_distance = 4, threshold = 0.1, voxelspacing = None, mask = slice(None)):
      r"""
      Takes a simple or multi-spectral image and returns the directionality of the image texture.
      It is just a value representing the strength of directionality, not the specific direction. 
      
      An edge detection is applied on the image. Then the edge strength and directional angle between
      the image axis are computed. A histogram of the directional angles is than used to calculate a
      qualitative value for directionality in ONE image layer. Note that there are n choose 2 layers
      in a n dimensional image. 
      
      Warning
      -------
      Experimental. There are still issues with finding the right maxs and mins in histogram and
      predefining the number of bins for the histogram.
  
      Parameters
      ----------
      image : array_like or list/tuple of array_like 
          A single image or a list/tuple of images (for multi-spectral case).
      voxelspacing : sequence of floats
          The side-length of each voxel.
      mask : array_like
          A binary mask for the image or a slice object
      min_distance : int
          minimal Distance between 2 local minima or maxima in the histogram. Default is 4.
      threshold : float
          Defines a threshold between 0 and 1. It is used to ignore angles of low edge strength
          in the histogram. Default is 0.1.
  
      Returns
      -------
      directionality : array
          Fdir is a value between 0 and 1. 1 represents a high directionality.
          Returns the directionality of an image in relation to one special image layer.
          The returned values are sorted like this. The axis are named v,w,x,y,z
          for a five dimensional image:
                                      w   x   y   z   v     x   y   z   v   w
          arctan(delta)| delta =    ---,---,---,---,---,  ---,---,---,---,---
                                      v   w   x   y   z     v   w   x   y   z
          There are always n choose k axis relations; n=image.ndim, k=2 (2 axis in every image layer).
          
      
      See Also
      --------
      
      """
      image = numpy.asarray(image)
      ndim = image.ndim
      # set default mask or apply given mask
      if not type(mask) is slice:
          if not type(mask[0] is slice):
              mask = numpy.array(mask, copy=False, dtype = numpy.bool)
      image = image[mask]
             
      # set default voxel spacing if not suppliec
      if None == voxelspacing:
          voxelspacing = tuple([1.] * ndim)
          
      if len(voxelspacing) != ndim:
          print("Voxel spacing and image dimensions do not fit.")
          return None
     
     # Calculate amount of combinations: n choose k, normalizing factor r and voxel spacing.    
      n = (factorial(ndim)/(2*factorial(ndim-2)))
      pi1_2 = numpy.pi/2.0
      r=1.0 / (pi1_2**2)
      vs = [slice(None,None,numpy.rint(ii)) for ii in voxelspacing]
     
      # Allocate memory, define constants
      Fdir = numpy.empty(n)
  
      # calculate differences by using Sobel-filter. (Maybe other filter kernel like Prewitt will do a better job)
      E = [sobel(image, axis=ndim-1-i) for i in range(ndim)]
      
      # The edge strength e(x,y) is used for thresholding.
      e = sum(E) / float(ndim)
      border = [numpy.percentile(e, 1),numpy.percentile(e, 99)]
      e[e < border[0]] = 0
      e[e > border[1]] = border[1]
      e -= border[0]
      e /= border[1]
      em = e > threshold
          
      for i in range(n):
          A = numpy.arctan((E[(i + (ndim+i)/ndim) % ndim][vs]) / (E[i%ndim][vs]+numpy.spacing(1))) # [0 , pi/2]
          A = A[em[vs]]
          # Calculate number of bins for the histogram. Watch out, this is just a work around! 
          # @TODO: Write a more stable code to prevent for minimum and maximum repetition when the same value in the Histogram appears multiple times in a row. Example: image = numpy.zeros([10,10]), image[:,::3] = 1
          bins = numpy.unique(A).size + min_distance        
          H = numpy.histogram(A, bins = bins, density=True)[0] # [0 , 1]
          H[H < numpy.percentile(H,1)] = 0.0
          H_peaks, H_valleys, H_range = find_valley_range(H)
          summe = 0.0
          for idx_ap in range(len(H_peaks)):
              for range_idx in range( H_valleys[idx_ap], H_valleys[idx_ap]+H_range[idx_ap]):
                  a=range_idx % len(H)
                  summe += (((pi1_2*a)/bins - (pi1_2 * H_peaks[idx_ap])/bins) **2) * H[a]
          Fdir[i] = 1.0 - r * summe 
          
      return Fdir