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:

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