docs: def cloudInCells(x, y, bins, weights=None):
      """
      Use cloud-in-cells binning algorithm. Only valid for equal-spaced linear bins.
      http://ta.twi.tudelft.nl/dv/users/Lemmens/MThesis.TTH/chapter4.html#tth_sEc2
      http://www.gnu.org/software/archimedes/manual/html/node29.html
  
      INPUTS:
          x: array of x-values
          y: array or y-values
          bins: [bins_x, bins_y] format, where bins_x corresponds to the bin edges along x-axis
          weights[None]: optionally assign a weight to each entry
      OUTPUTS:
          histogram:
          bins_x:
          bins_y:
      """
  
      # For consistency, the variable names should be changed in this function, but low priority...
      
      x_bins = np.array(bins[0])
      delta_x = x_bins[1] - x_bins[0]
      # Overflow and underflow bins
      x_bins = np.insert(x_bins, 0, x_bins[0] - delta_x)
      x_bins = np.append(x_bins, x_bins[-1] + delta_x)
      y_bins = np.array(bins[1])
      delta_y = y_bins[1] - y_bins[0]
      y_bins = np.insert(y_bins, 0, y_bins[0] - delta_y)
      y_bins = np.append(y_bins, y_bins[-1] + delta_y)
  
      x_bound_cut = np.logical_and(x >= x_bins[0], x <= x_bins[-1])
      y_bound_cut = np.logical_and(y >= y_bins[0], y <= y_bins[-1])
      bound_cut = np.logical_and(x_bound_cut, y_bound_cut)
  
      if not np.any(weights):
          bound_weights = np.ones(len(x))[bound_cut]
      else:
          bound_weights = np.array(weights)[bound_cut]
  
      x_vals = np.array(x)[bound_cut]
      y_vals = np.array(y)[bound_cut]
  
      x_width = x_bins[1] - x_bins[0]
      y_width = y_bins[1] - y_bins[0]
  
      x_centers = x_bins[0: -1] + (0.5 * x_width)  
      y_centers = y_bins[0: -1] + (0.5 * y_width)
  
      dx = x_vals - x_centers[np.digitize(x_vals, x_bins) - 1]
      dy = y_vals - y_centers[np.digitize(y_vals, y_bins) - 1]
  
      ux = ((dx / x_width) * (dx >= 0)) +\
           ((1. + (dx / x_width)) * (dx < 0))
      lx = 1. - ux
  
      uy = ((dy / y_width) * (dy >= 0)) +\
           ((1. + (dy / y_width)) * (dy < 0))
      ly = 1. - uy
  
      new_x_vals = []
      new_y_vals = []
      cell_weights = []
  
      # 4 corners
      new_x_vals.append(x_vals + (0.5 * x_width))
      new_y_vals.append(y_vals + (0.5 * y_width))
      cell_weights.append(bound_weights * ux * uy)
  
      new_x_vals.append(x_vals + (0.5 * x_width))
      new_y_vals.append(y_vals - (0.5 * y_width))
      cell_weights.append(bound_weights * ux * ly)
      
      new_x_vals.append(x_vals - (0.5 * x_width))
      new_y_vals.append(y_vals + (0.5 * y_width))
      cell_weights.append(bound_weights * lx * uy)
      
      new_x_vals.append(x_vals - (0.5 * x_width))
      new_y_vals.append(y_vals - (0.5 * y_width))
      cell_weights.append(bound_weights * lx * ly)
  
      new_x_vals = np.concatenate(new_x_vals)
      new_y_vals = np.concatenate(new_y_vals)
      cell_weights = np.concatenate(cell_weights)
  
      result = np.histogram2d(new_x_vals, new_y_vals,
                                 bins = [x_bins, y_bins],
                                 weights = cell_weights)[0]
   
      result = np.transpose(result[1: result.shape[0] - 1])[1: result.shape[1] - 1]
  
      return result, x_bins, y_bins