Efficient point grouping algorithm
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})&lt;\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)&lt;\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.radii  = radius
        self.nodelst = nodelst #nodelst is array of coordinate tuples, graph contains indices as nodes
        self.__make_edges__()

    def __make_edges__(self):
        self.igraph.add_edges_from( sp.cKDTree(self.nodelst).query_pairs(r=self.radii) )

    def get_conn_comp(self):
        ind = [list(x) for x in nx.connected_components(self.igraph) if len(x)&gt;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).