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