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