def coarseness(image, voxelspacing = None, mask = slice(None)):
    r"""
    Takes a simple or multi-spectral image and returns the coarseness of the texture.
    
    Step1  At each pixel, compute six averages for the windows of size 2**k x 2**k,
            k=0,1,...,5, around the pixel. 
    Step2  At each pixel, compute absolute differences E between the pairs of non 
            overlapping averages in every directions.
    step3  At each pixel, find the value of k that maximises the difference Ek in either 
            direction and set the best size Sbest=2**k
    step4  Compute the coarseness feature Fcrs by averaging Sbest over the entire image.

    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
        
    Returns
    -------
    coarseness : float
        The size of coarseness of the given texture. It is basically the size of
        repeating elements in the image. 
        
    See Also
    --------
    
    
    """
    # Step1:  At each pixel (x,y), compute six averages for the windows
    # of size 2**k x 2**k, k=0,1,...,5, around the pixel.

    image = numpy.asarray(image, dtype=numpy.float32)
   
  
    # 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.] * image.ndim)
    
    if len(voxelspacing) != image.ndim:
        print("Voxel spacing and image dimensions do not fit.")
        return None
    # set padding for image border control
    padSize = numpy.asarray([(numpy.rint((2**5.0) * voxelspacing[jj]),0) for jj in range(image.ndim)]).astype(numpy.int)
    Apad = numpy.pad(image,pad_width=padSize, mode='reflect')

    # Allocate memory
    E = numpy.empty((6,image.ndim)+image.shape)

    # prepare some slicer 
    rawSlicer           = [slice(None)] * image.ndim
    slicerForImageInPad = [slice(padSize[d][0],None)for d in range(image.ndim)]

    for k in range(6):

        size_vs = tuple(numpy.rint((2**k) * voxelspacing[jj]) for jj in range(image.ndim))
        A = uniform_filter(Apad, size = size_vs, mode = 'mirror')

        # Step2: At each pixel, compute absolute differences E(x,y) between 
        # the pairs of non overlapping averages in the horizontal and vertical directions.
        for d in range(image.ndim):
            borders = numpy.rint((2**k) * voxelspacing[d])
            
            slicerPad_k_d   = slicerForImageInPad[:]
            slicerPad_k_d[d]= slice((padSize[d][0]-borders if borders < padSize[d][0] else 0),None)
            A_k_d           = A[slicerPad_k_d]

            AslicerL        = rawSlicer[:]
            AslicerL[d]     = slice(0, -borders)
            
            AslicerR        = rawSlicer[:]
            AslicerR[d]     = slice(borders, None)

            E[k,d,...] = numpy.abs(A_k_d[AslicerL] - A_k_d[AslicerR])

    # step3: At each pixel, find the value of k that maximises the difference Ek(x,y)
    # in either direction and set the best size Sbest(x,y)=2**k
    
    k_max = E.max(1).argmax(0)
    dim = E.argmax(1)
    dim_vox_space = numpy.asarray([voxelspacing[dim[k_max.flat[i]].flat[i]] for i in range(k_max.size)]).reshape(k_max.shape) 
    S = (2**k_max) * dim_vox_space

    # step4: Compute the coarseness feature Fcrs by averaging Sbest(x,y) over the entire image.
    return S.mean()