# bin 2d array such that each bin contains equal number of samples?

I'm trying to bin a 2d array of points such that each bin has an equal number of samples contained in it. This would mean having bins of differing ranges and possibly shapes.

My question is: can scipy.spatial.cKDTree somehow help solve this problem? It doesn't seem clear to me how that could be achieved.

You could compute the Voronoi diagram. Each cell in this diagram contains exactly one point and taken together the cells cover the domain. If you need, say, $$n$$ points in each bin, it is then a matter of picking a cell, merging it with its $$n-1$$ neighbors, and doing the same for the remaining cells.

• How do you deal with a cell that already is the neighbour of another cell? When I tried implementing a type of nearest neighbour approach, that involved creating spaced points and then finding the n nearest neighbours from those points to points in the dataset that weren't already used. It didn't give a desirable result as the last remaining cells were forced to take points that were very spread out. I can't think of a way to avoid that when using a Voronoi diagram. Feb 16 '21 at 22:14
• You want the bins to be connected? I’m not sure that is always possible. What is the property, apart from having the same number of points, that you are looking for in the bins? Feb 17 '21 at 9:48
• I suppose I wasn't very clear in my original question, apologies. I have a vague idea that points of the same bin ought to be close to other points of the same bin. One point entirely encapsulated by points of a different bin, would be an undesirable result. Ideally, every bin could be described by a simple polygon that isn't unreasonably elongated. Bins of vastly different sizes (due to large differences in densities in the data) are fine as long as the number of points in them remains constant. Thank you for your patience. Feb 17 '21 at 12:42
• Then using a kdtree sounds to me like the way to go. I’m afraid I don’t know the function in scipy. Feb 18 '21 at 8:14
• Ok, thank you! One more question: is a "balanced" kdtree needed for the algorithm that I described? Feb 18 '21 at 10:20

Using the .tree attribute of cKDTree gives you direct access to the kdtree. You can recursively navigate it like this:

from scipy.spatial import cKDTree
tree = cKDTree(points) # let points be an array of shape (n,2)
groups = []
depth = 5 # depth of the tree, ie: depth=5 results in 64 bins
def recurse(node, i=0):
g = node.greater
l = node.lesser
if(i == depth):
groups.append(g.data_points)
groups.append(l.data_points)
return
recurse(g, i=i+1)
recurse(l, i=i+1)

recurse(tree.tree)


The variable "groups" will now be a list of lists where groups[a][b] denotes the ath bin and bth point inside that bin. These bins will most of the time contain the same number of points (within 1).

There is an issue that I believe has to do with the fact the a KDTree only splits groups perpendicular to the coordinate axes. This means that if your data has a grid-like structure (that is to say there are many points that share an x or y coordinate), it is possible that groups will have different sizes. While there might be a better solution to this problem, I don't know what it is. So I present to you my very hacky solution:

Rotate the grid by some small angle (10°), then compute the tree, then rotate the partitions back. In hopes that the data doesn't have some kind of symmetry at 10°. Please make a suggestion if you have a better idea, here is a full example of some sample code:

import numpy as np
from matplotlib import pyplot as plt
from scipy.spatial import cKDTree

# Generate example data
h = 30  # height
w = 30  # width
l = 3  # layers
n_points = l * h * w
points = np.zeros((n_points, 2))
for i in range(n_points):
hx = i % h
wy = int(i/w) % w
points[i] = [1+hx + np.random.normal(scale=0.1),
1 + wy + np.random.normal(scale=0.1)]

# Setup rotation matrix
deg = 30
j = np.matrix([[c, s], [-s, c]])

# Generate KDTree with points rotated
tree = cKDTree((j @ points.T).T,  balanced_tree=True)

# Navigate KDTree
groups = []
def recurse(node, i=0):
g = node.greater
l = node.lesser
if(i == 5):
groups.append(g.data_points)
groups.append(l.data_points)
return
recurse(g, i=i+1)
recurse(l, i=i+1)

recurse(tree.tree)

# Setup rotation matrix for -10°
j = np.matrix([[c, s], [-s, c]])

# Count the size of bins, rotate points back and plot bins.
list_len = []
for group in groups:
points_rot = (j @ group.T).T
list_len.append(len(group))
plt.plot(points_rot[:, 0], points_rot[:, 1], 'o', markersize=2)
print(f"There are {len(groups)} groups.", end=' ')
print(f"biggest - smallest = {max(list_len)} - {min(list_len)} = {max(list_len) - min(list_len)}")

plt.show()


There are 64 groups. biggest - smallest = 43 - 42 = 1