docs: def sum_of_entropies(image, mask, binary_image):
      """Bin the foreground and background pixels and compute the entropy 
      of the distribution of points among the bins
      """
      mask=mask.copy()
      mask[np.isnan(image)] = False
      if not np.any(mask):
          return 0
      #
      # Clamp the dynamic range of the foreground
      #
      minval = np.max(image[mask])/256
      if minval == 0:
          return 0
      clamped_image = image.copy()
      clamped_image[clamped_image < minval] = minval
      #
      # Smooth image with -8 bits of noise
      #
      image = smooth_with_noise(clamped_image, 8)
      im_min = np.min(image)
      im_max = np.max(image)
      #
      # Figure out the bounds for the histogram
      #
      upper = np.log2(im_max)
      lower = np.log2(im_min)
      if upper == lower:
          # All values are the same, answer is log2 of # of pixels
          return math.log(np.sum(mask),2) 
      #
      # Create log-transformed lists of points in the foreground and background
      # 
      fg = image[binary_image & mask]
      bg = image[(~ binary_image) & mask]
      if len(fg) == 0 or len(bg) == 0:
          return 0
      log_fg = np.log2(fg)
      log_bg = np.log2(bg)
      #
      # Make these into histograms
      hfg = numpy_histogram(log_fg, 256, range=(lower,upper))[0]
      hbg = numpy_histogram(log_bg, 256, range=(lower,upper))[0]
      #hfg = scipy.ndimage.histogram(log_fg,lower,upper,256)
      #hbg = scipy.ndimage.histogram(log_bg,lower,upper,256)
      #
      # Drop empty bins
      #
      hfg = hfg[hfg>0]
      hbg = hbg[hbg>0]
      if np.product(hfg.shape) == 0:
          hfg = np.ones((1,),int)
      if np.product(hbg.shape) == 0:
          hbg = np.ones((1,),int)
      #
      # Normalize
      #
      hfg = hfg.astype(float) / float(np.sum(hfg))
      hbg = hbg.astype(float) / float(np.sum(hbg))
      #
      # Compute sum of entropies
      #
      return np.sum(hfg * np.log2(hfg)) + np.sum(hbg*np.log2(hbg))