Efficient point grouping algorithm

I'm trying to find a good solution to the following problem: Given a set of 3d real number coordinates inside a cube, I want to group them into containers. The containers should contain all elements that are in predefined distance to at least one of the other elements. Should two elements be closer than a predefined distance, they are to end in the same container. Note: Two elements that are further apart may end in the same container e.g. if they are both close to a different element

The most straight forward solution I can think of is to go through all points, and compare them with all other points that haven't been visited before to find out if they are in proximity, but it seems rather inefficient.

I tried to find a solution online, but (maybe because I'm missing the correct keywords) can't find any. I stumbled upon k-d-trees and similar concepts but as far as I understood they don't help much here. Are there good algorithms out there that can take care of this or does anybody have a suggestion on what I could do?

Below is a (badly hand drawn) illustration of a possible 2d case with reference distance d marked in blue and the two groups I want to end up with in red. Thanks in advance for any help.

Edit: as pointed out in the comments, this problem can be thought of as finding connected components of a graph. The graph nodes are the points from the data set and two nodes are adjacent if their distance is smaller than the reference value. The problem I see is that establishing the edges seems to be O(n^2), because I don't see a way that doesn't require me to do pairwise checks on all points of the set, is that correct?

• How many groups? 2 or more ? Dec 26, 2017 at 15:56
• Unforunately unable to know that beforehand :/ The requirement is also strict. Meaning that, if there is a point that is closer than distance d to some other point, they have to end up in the same group. Edit: The distance is fixed at all times. I don't want to vary the distance to end up with better clustering or something like that. Dec 26, 2017 at 16:01
• Still it’s vague. What if you have points at $(1,0,0)$, $(2,0,0)$ and $(3,0,0)$. Your $d=1.1$ and two groups. What’s the solution? Dec 26, 2017 at 17:14
• If d is 1.1 and I call the points you mentioned A,B,C, then A is in proximity of B, B in proximity of C, so they all should end up in one and the same class {A,B,C}. I dont specify the number of groups beforehand. A fourth point (5,0,0) is close to neither of A,B,C and should thus be in a new class of his own. Dec 27, 2017 at 2:15
• sounds exactly like connected components to me Dec 27, 2017 at 6:06

There is an efficient algorithm whose running time is nearly linear in most cases.

Suppose you the distance threshold is $d$, i.e., you want two elements at distance $\le d$ to be in the same group.

Then you can divide up space into cubes of side length $2d$, and store each element in the cube it is contained in. The cubes can be stored in a hash table for efficient lookup. For instance, given a point at coordinates $(x,y,z)$, it is associated with the cube whose lower-left corner is at coordinates $(\lfloor x/2d \rfloor, \lfloor y/2d \rfloor, \lfloor z/2d \rfloor \rfloor)$; you can hash those coordinates and store the point in the corresponding bucket of the hashtable.

Now, if two points are at distance $\le d$ apart, then they must be in adjacent cubes. So, you can find all points that need to be put into the same group by looking at pairs of adjacent cubes. In particular, the algorithm becomes:

• For each cube $C$:
• For each cube $C'$ that is adjacent to $C$:
• For each point $P$ in $C$ and each point $P'$ in $C'$:
• If the distance between $P,P'$ is at most $d$, merge them into the same group.

You can keep track of the groups using a Union-Find data structure. With this approach, the running time will be $O(n \alpha(n))$, assuming that each cube has at most a constant number of points in it (which will typically be the case for most arrangements of points). Here we've used the fact that each cube is adjacent to only 26 other cubes, so each point $P$ only needs to be compared to a constant number of other points.

• Thank you for your answer! I stumbled upon this algorithm after someone mentioned you could build the union disk graph in linear time in this paper: sciencedirect.com/science/article/pii/0020019077900709 I tried to implement it, and it works, but DBSCAN runs faster on my machine. I assume it is because my implementation is not really optimal. Do you happen to know how the algorithm supposedly compares to DBSCAN? Dec 28, 2017 at 5:12
• @Banana, I don't, but perhaps you'd like to write an answer mentioning those two algorithms, and giving a full citation to the paper (title, where published, authors)? Sounds useful to me. Does DBSCAN really solve your problem? Doesn't it use a different set of criteria for grouping elements? Also, I suspect the implementation uses k-d trees or similar, which you rejected (though to be honest I'm not clear on why you rejected them).
– D.W.
Dec 28, 2017 at 18:24
• I will write a little about it later :) I didn't really reject k-d-trees, it was just that I don't properly understand them (thats why I did not see how I could utilize them) and how to implement them. Actually they kind of appear in my final implementation. Dec 29, 2017 at 4:46
• added a lengthy answer, hope it is indeed useful for you. Dec 29, 2017 at 8:31

My Solutions

I will try to elaborate on what solutions I found, their implementations and those implementation's running time. As I don't have any real computer science background, please point out any mistakes/unclarities you find.

First: The problem, in general terms, is to partition a set of points from $A \subset \mathbb{R}^d$ into groups such that every point in them is $\epsilon$-reachable from every other point in the group. Reachable,for $x_1, x_n \in A$, means there is a sequence of points $x_1,\dots,x_n \in A$ s.t. $d(x_i,x_{i+1})<\epsilon$. In other terms: If two points are closer than $\epsilon$ to each other, the algorithm must assign them to the same group/cluster.

DBSCAN: DBSCAN is a super general algorithm. I will use $B_{\epsilon}(x)$ for the $d$-dimensional sphere around $x$ with respect to some arbitrary metric $d(\cdot,\cdot)$. DBSCAN finds clusters by defining:

• $x \in A$ is a core point, if at least $N$ other points of $A$ lie in $B_{\epsilon}(x)$
• $x \in A$ is directly reachable from $y \in A$ if both are core points and $d(x,y)<\epsilon$
• $x \in A$ is reachable from $y \in A$ if there exists a sequence $x_1 =x, \dots, x_n =y \in A$ s.t. $x_{i+1}$ is reachable from $x_i$
• non reachable points are outliers

Now, as you might see, for $N=2$ (if you include every point itself when counting) and $d(\cdot, \cdot)$ the usual euclidean metric this defines exactly the problem I stated above. DBSCAN clusters together core-points deterministically, but also adds non-core points to clusters. So if $x \in A$ is reachable from core point $y$, it might end up in the same cluster without being a core point itself. I say might, because non-core points may be reachable from different clusters. The assignation is then determined by the order in which the algorithm runs through the clusters. This ambiguity however resolves with $N=2$, as then every point that is reachable from another is automatically a core point, that is reachability as stated above is then a symmetric relation. Clusters are then unambiguous. In case of an indexing structure that does neighborhood querys in $O(logn)$, the runtime complexity is apparently $O(n \; logn)$, but don't take my word for it.

I used scikit-learn's implementation of DBSCAN (http://scikit-learn.org/stable/modules/generated/sklearn.cluster.DBSCAN.html). I failed to find out what runtime complexity this implementation has, but the docs state memory complexity of $O(n.d)$.

Intersection graph: As it has been proposed in a different answer: One can also think of the problem as a graph-theoretical one. The graphs vertices are represented by the points in $A$, with two vertices being adjacent iff their points are in $\epsilon$-distance of each other. On this graph the problem translates to finding the connected components. The paper of Bentley, Jon L.; Stanat, Donald F.; Williams, E. Hollins, Jr. (1977), "The complexity of finding fixed-radius near neighbors", Information Processing Letters, 6 (6): 209–212 provides an algorithm for constructing this graph that runs in linear time and is essentially identical to the one proposed by @D.W.. Finding the connected components can also be done in linear time, e.g. by breadth- or depth-first algortihms.

Python Code: For the Code, I mainly used the DBSCAN implementation of scikit-learn, networkx for graphs and finding connected components (it is open source but I couldnt find any reference on the runtime complexity) and scipy.spatial's cKDTree implementation (which I do not understand in detail).

from sklearn.cluster import DBSCAN
from sklearn.datasets.samples_generator import make_blobs
import networkx as nx
import scipy.spatial as sp

def cluster(data, epsilon,N): #DBSCAN, euclidean distance
db     = DBSCAN(eps=epsilon, min_samples=N).fit(data)
labels = db.labels_ #labels of the found clusters
n_clusters = len(set(labels)) - (1 if -1 in labels else 0) #number of clusters
clusters   = [data[labels == i] for i in range(n_clusters)] #list of clusters
return clusters, n_clusters


Is my minimal implementation of the previously described goal. The method returns clusters and their size. For the graph method, I wrote a little class. The class describes a graph whose vertices are the indices of the corresponding points in $A$ as list. Using the cKDTree's data structure, finding pairs of reachable points is lightning fast. You can then slice out the elements via the indices.

class IGraph:
def __init__(self, nodelst=[], radius = 1):
self.igraph = nx.Graph()
self.nodelst = nodelst #nodelst is array of coordinate tuples, graph contains indices as nodes
self.__make_edges__()

def __make_edges__(self):

def get_conn_comp(self):
ind = [list(x) for x in nx.connected_components(self.igraph) if len(x)>1]
return [self.nodelst[indlist] for indlist in ind]

def graph_cluster(data, epsilon):
graph = IGraph(nodelst = data, radius = epsilon)
clusters = graph.get_conn_comp()
return clusters, len(clusters)


I also did a little benchmark, but it is in no way representative. I create a cluster with:

centers = [[1, 1,1], [-1, -1,1], [1, -1,1]]
X,_ = make_blobs(n_samples=N, centers=centers, cluster_std=0.4,
random_state=0)


For different N, the methods compare differently, on my machine for example:

• $N = 20000$: "cluster": 293.711 ms, "graph_cluster": 422.089 ms
• $N = 7000$: "cluster": 87.455 ms, "graph_cluster": 95.322 ms
• $N = 700$: "cluster": 3.764 ms, "graph_cluster": 1.689 ms

For small datasets, and this specific clustering, the graph solution wins. But apparently the DBSCAN shows better scaling behavior. I haven't looked at different dimensions of the data, but that may also be interesting. Also: the results are exactly the same for both algortihms. Additional Note: The graph solution scales heavily with $\epsilon$, which is 0.1 in my code. Larger values make it far worse! In other words: If you have small clusters, the graph solution outperforms DBSCAN by ages, for large clusters it is the other way around (apparently).

You can build the graph according to the following rules:

1) nodes corresponds to points

2) two nodes are connected if and only if the distance between points is less than $d$

Then you either look for connected components in the graph (classical and very fast, use DFS) or apply another graph clustering/partitioning algorithm. If splitting into two groups MIN-CUT is a reasonable choice.

• Yes that is exactly what I thought. But identifying the graph first and then running one of those algorithms should have a time complexity of worse than O(n^2), right? Dec 27, 2017 at 2:19
• but its a good start thank you, I will take a look at the algorithms Dec 27, 2017 at 2:20
• @banana also, I might have misunderstood. You can construct the unit disk graph with given centers in linear time, if I remember correctly. Then DFS should be extremely fast ($O(m+n)$ though). Please check if a unit disk graph is a valid model Dec 27, 2017 at 6:40
• I didn't know you could construct the unit disk graph in linear time, that sounds perfect! I will check that out and see how it works out. Also I stumbled upon the DBSCAN which can take on the task (epsilon = d and minPts = 2) and seems to promose O(n*log(n)) Dec 27, 2017 at 6:45