docs: def entropy(data, bins=256):
      """Compute a threshold using Ray's entropy measurement
      
      data           - an array of intensity values between zero and one
      bins           - we bin the data into this many equally-spaced bins, then pick
                       the bin index that optimizes the metric
      """
      
      data = np.atleast_1d(data)
      data = data[~ np.isnan(data)]
      if len(data) == 0:
          return 0
      elif len(data) == 1:
          return data[0]
  
      if bins > len(data):
          bins = len(data)
      data.sort()
      var = running_variance(data)+1.0/512.0
      rvar = np.flipud(running_variance(np.flipud(data)))+1.0/512.0 
      thresholds = data[1:len(data):len(data)//bins]
      w = np.arange(0,len(data)-1,len(data)//bins)
      score_low = w * np.log(var[0:len(data)-1:len(data)//bins] *
                             w * np.sqrt(2*np.pi*np.exp(1)))
      score_low[np.isnan(score_low)]=0
      
      w = len(data) - np.arange(1,len(data),len(data)//bins)
      score_high = w * np.log(rvar[1:len(data):len(data)//bins] * w *
                              np.sqrt(2*np.pi*np.exp(1)))
      score_high[np.isnan(score_high)]=0
      scores = score_low + score_high
      index = np.argwhere(scores == scores.min()).flatten()
      if len(index)==0:
          return 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 (thresholds[index_low]+thresholds[index_high]) / 2