docs: def backgr(img, mask = None, mode=MODE_AUTO, thresh=2, splinepoints=5, scale=1, 
             maxiter=40, convergence = .001):
      '''Iterative spline-based background correction.
      
      mode - one of MODE_AUTO, MODE_DARK, MODE_BRIGHT or MODE_GRAY
      thresh - thresh is threshold to cut at, in units of sigma. 
      splinepoints - # of points in spline in each direction
      scale - scale the image by this factor (e.g. 2 = operate on 1/2 of the points)
      maxiter - maximum # of iterations
      convergence - result has converged when the standard deviation of the
                    difference between iterations is this fraction of the
                    maximum image intensity.
  Spline mesh is splinepoints x splinepoints. Modes are
  defined by the background intensity, Larger thresh->slower but
  more stable convergence. Returns background matrix.'''
  
      assert img.ndim == 2, "Image must be 2-d"
      assert splinepoints >= 3, "The minimum grid size is 3x3"
      assert maxiter >= 1
      assert mode in [MODE_AUTO, MODE_BRIGHT, MODE_DARK, MODE_GRAY], mode + " is not a valid background mode"
      orig_shape = np.array(img.shape).copy()
      input_mask = mask
      if mask is None:
          mask = np.ones(orig_shape, dtype=bool) #start with mask = whole image
          clip_imin = clip_jmin = 0
          clip_imax = img.shape[0]
          clip_jmax = img.shape[1]
          clip_shape = orig_shape
      else:
          isum = np.sum(mask, 1)
          jsum = np.sum(mask, 0)
          clip_imin = np.min(np.argwhere(isum != 0))
          clip_imax = np.max(np.argwhere(isum != 0)) + 1
          clip_jmin = np.min(np.argwhere(jsum != 0))
          clip_jmax = np.max(np.argwhere(jsum != 0)) + 1
          clip_shape = np.array([clip_imax - clip_imin, clip_jmax - clip_jmin])
          
      subsample_shape = (clip_shape / scale).astype(int)
      ratio = (clip_shape.astype(float) - 1) / (subsample_shape.astype(float) - 1)
      transform = np.array([[ratio[0], 0], [0, ratio[1]]])
      inverse_transform = np.array([[1.0 / ratio[0], 0],
                                    [0, 1.0 / ratio[1]]])
      
      img = affine_transform(img[clip_imin:clip_imax, clip_jmin:clip_jmax],
                             transform, output_shape=tuple(subsample_shape),
                             order = 2)
      mask = affine_transform(mask[clip_imin:clip_imax, 
                                   clip_jmin:clip_jmax].astype(float), 
                              transform, 
                              output_shape=tuple(subsample_shape),
                              order = 2) > .5
      orig_mask = mask
      
      if mode=='auto':
          mode = automode(img[orig_mask])
      elif mode=='dark' or mode=='low':
          mode = -1
      elif mode=='bright' or mode=='high':
          mode = 1
      elif mode=='gray' or mode=='grey' or mode=='mid':
          mode = 0
  
      # Base the stop criterion on a fraction of the image dynamic range
      
      stop_criterion=max((np.max(img) - np.min(img)) * convergence,
                         np.finfo(img.dtype).eps)
      [r,c] = img.shape
  
      oldres = np.zeros((r,c)) #old background
  
      for i in range(maxiter):
          px, py, pz = splineimage(img, splinepoints, np.array(mask)) #now with mask
          res = evalspline2d(np.arange(c), np.arange(r), px, py, pz)
          comp = img - res
  
          diff = res[orig_mask] - oldres[orig_mask]
          ###Compute std. deviation in same way as matlab std(), (not numpy.std())
          stddiff = unbiased_std(diff)
  
          if stddiff < stop_criterion:#stop_criterion instead of .004
              break
          elif i==maxiter:
              warnings.warn('Background did not converge after %d iterations.\nMake sure that the foreground/background mode is correct.'%(i))
  
          oldres = res
          
          #calculate new mask
          backgr = comp[mask]
          sigma = unbiased_std(backgr)
          cut = sigma * thresh
          
          if mode < 0:
              mask = comp < cut
          elif mode > 0:
              mask = comp > -cut
          else:
              mask = abs(comp) < cut
          mask &= orig_mask
          nnz = np.sum(mask)
          if nnz < .01 * np.sum(orig_mask):
              warnings.warn('Less than 1%% of the pixels used for fitting,\ntry starting again with a larger threshold value')
              break
      
      output = np.zeros(orig_shape, img.dtype)
      output[clip_imin:clip_imax, clip_jmin:clip_jmax] = \
            affine_transform(res, inverse_transform, 
                             output_shape = tuple(clip_shape),
                             order = 3)
      if input_mask is not None:
          output[~input_mask] = 0
      return output