docs: def otsu(data, min_threshold=None, max_threshold=None,bins=256):
      """Compute a threshold using Otsu's method
      
      data           - an array of intensity values between zero and one
      min_threshold  - only consider thresholds above this minimum value
      max_threshold  - only consider thresholds below this maximum value
      bins           - we bin the data into this many equally-spaced bins, then pick
                       the bin index that optimizes the metric
      """
      assert min_threshold is None or max_threshold is None or min_threshold < max_threshold
      def constrain(threshold):
          if not min_threshold is None and threshold < min_threshold:
              threshold = min_threshold
          if not max_threshold is None and threshold > max_threshold:
              threshold = max_threshold
          return threshold
      
      data = np.atleast_1d(data)
      data = data[~ np.isnan(data)]
      if len(data) == 0:
          return (min_threshold if not min_threshold is None
                  else max_threshold if not max_threshold is None
                  else 0)
      elif len(data) == 1:
          return constrain(data[0])
      if bins > len(data):
          bins = len(data)
      data.sort()
      var = running_variance(data)
      rvar = np.flipud(running_variance(np.flipud(data))) 
      thresholds = data[1:len(data):len(data)//bins]
      score_low = (var[0:len(data)-1:len(data)//bins] * 
                   np.arange(0,len(data)-1,len(data)//bins))
      score_high = (rvar[1:len(data):len(data)//bins] *
                    (len(data) - np.arange(1,len(data),len(data)//bins)))
      scores = score_low + score_high
      if len(scores) == 0:
          return constrain(thresholds[0])
      index = np.argwhere(scores == scores.min()).flatten()
      if len(index)==0:
          return constrain(thresholds[0])
      #
      # Take the average of the thresholds to either side of
      # the chosen value to get an intermediate in cases where there is
      # a steep step between the background and foreground
      index = index[0]
      if index == 0:
          index_low = 0
      else:
          index_low = index-1
      if index == len(thresholds)-1:
          index_high = len(thresholds)-1
      else:
          index_high = index+1 
      return constrain((thresholds[index_low]+thresholds[index_high]) / 2)