diff --git a/.travis.yml b/.travis.yml
index 6dac1059..d2d62fff 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -25,9 +25,10 @@ addons:
     - ubuntu-toolchain-r-test
     - sourceline: 'deb http://downloads.skewed.de/apt/trusty trusty universe'
     packages:
+    - libqt5gui5  # pyqt5>5.11 otherwise cannot load the xcb platform plugin
+    - libflann-dev
     - python3-graph-tool
     - python-graph-tool
-    - libqt5gui5  # pyqt5>5.11 fails to load the xcb platform plugin without it
 
 install:
   - pip install --upgrade --upgrade-strategy eager .[dev]
diff --git a/README.rst b/README.rst
index 12f8e899..8e1c307d 100644
--- a/README.rst
+++ b/README.rst
@@ -39,6 +39,15 @@ A (mostly unmaintained) `Matlab version <https://lts2.epfl.ch/gsp>`_ exists.
 .. |conda| image:: https://anaconda.org/conda-forge/pygsp/badges/installer/conda.svg
    :target: https://anaconda.org/conda-forge/pygsp
 
+
+The PyGSP is a Python package to ease
+`Signal Processing on Graphs <https://arxiv.org/abs/1211.0053>`_.
+The documentation is available on
+`Read the Docs <https://pygsp.readthedocs.io>`_
+and development takes place on
+`GitHub <https://github.com/epfl-lts2/pygsp>`_.
+A (mostly unmaintained) `Matlab version <https://epfl-lts2.github.io/gspbox-html>`_ exists.
+
 The PyGSP facilitates a wide variety of operations on graphs, like computing
 their Fourier basis, filtering or interpolating signals, plotting graphs,
 signals, and filters. Its core is spectral graph theory, and many of the
diff --git a/doc/history.rst b/doc/history.rst
index 11b5286c..77be17cc 100644
--- a/doc/history.rst
+++ b/doc/history.rst
@@ -24,6 +24,9 @@ History
 * New implementation of the Sensor graph that is simpler and scales better.
 * A new learning module with three functions to solve standard semi-supervised
   classification and regression problems.
+* A much improved, fixed, documented, and tested NNGraph. The user can now
+  select the backend and similarity kernel. The radius can be estimated and
+  features standardized. (PR #43)
 * Import and export graphs and their signals to NetworkX and graph-tool.
 * Save and load graphs and theirs signals to / from GraphML, GML, and GEXF.
 * Documentation: path graph linked to DCT, ring graph linked to DFT.
@@ -31,6 +34,7 @@ History
   taste of what the library can do, and to start working from a code snippet.
 * Merged all the extra requirements in a single dev requirement.
 
+
 Experimental filter API (to be tested and validated):
 
 * evaluate a filter bank with ``g(values)``
diff --git a/doc/tutorials/optimization.rst b/doc/tutorials/optimization.rst
index 7127a9c5..b3740d85 100644
--- a/doc/tutorials/optimization.rst
+++ b/doc/tutorials/optimization.rst
@@ -85,7 +85,7 @@ We start with the graph TV regularization. We will use the :class:`pyunlocbox.so
     >>> prob1 = pyunlocbox.solvers.solve([d, r, f], solver=solver,
     ...                                  x0=x0, rtol=0, maxit=1000)
     Solution found after 1000 iterations:
-        objective function f(sol) = 2.250584e+02
+        objective function f(sol) = 2.256055e+02
         stopping criterion: MAXIT
     >>>
     >>> fig, ax = G.plot(prob1['sol'])
@@ -107,7 +107,7 @@ This figure shows the label signal recovered by graph total variation regulariza
     >>> prob2 = pyunlocbox.solvers.solve([r, f], solver=solver,
     ...                                  x0=x0, rtol=0, maxit=1000)
     Solution found after 1000 iterations:
-        objective function f(sol) = 6.504290e+01
+        objective function f(sol) = 4.376481e+01
         stopping criterion: MAXIT
     >>>
     >>> fig, ax = G.plot(prob2['sol'])
diff --git a/pygsp/_nearest_neighbor.py b/pygsp/_nearest_neighbor.py
new file mode 100644
index 00000000..15dd97ce
--- /dev/null
+++ b/pygsp/_nearest_neighbor.py
@@ -0,0 +1,247 @@
+# -*- coding: utf-8 -*-
+
+from __future__ import division
+
+import numpy as np
+from scipy import sparse, spatial
+from pygsp import utils
+
+_logger = utils.build_logger(__name__)
+
+def _scipy_pdist(features, metric, order, kind, k, radius, params):
+    if params:
+        raise ValueError('unexpected parameters {}'.format(params))
+    metric = 'cityblock' if metric == 'manhattan' else metric
+    metric = 'chebyshev' if metric == 'max_dist' else metric
+    params = dict(metric=metric)
+    if metric == 'minkowski':
+        params['p'] = order
+    dist = spatial.distance.pdist(features, **params)
+    dist = spatial.distance.squareform(dist)
+    if kind == 'knn':
+        neighbors = np.argsort(dist)[:, :k+1]
+        distances = np.take_along_axis(dist, neighbors, axis=-1)
+    elif kind == 'radius':
+        distances = []
+        neighbors = []
+        for distance in dist:
+            neighbor = np.flatnonzero(distance < radius)
+            neighbors.append(neighbor)
+            distances.append(distance[neighbor])
+    return neighbors, distances
+
+
+def _scipy_kdtree(features, _, order, kind, k, radius, params):
+    if order is None:
+        raise ValueError('invalid metric for scipy-kdtree')
+    eps = params.pop('eps', 0)
+    tree = spatial.KDTree(features, **params)
+    params = dict(p=order, eps=eps)
+    if kind == 'knn':
+        params['k'] = k + 1
+    elif kind == 'radius':
+        params['k'] = None
+        params['distance_upper_bound'] = radius
+    distances, neighbors = tree.query(features, **params)
+    return neighbors, distances
+
+
+def _scipy_ckdtree(features, _, order, kind, k, radius, params):
+    if order is None:
+        raise ValueError('invalid metric for scipy-kdtree')
+    eps = params.pop('eps', 0)
+    tree = spatial.cKDTree(features, **params)
+    params = dict(p=order, eps=eps, n_jobs=-1)
+    if kind == 'knn':
+        params['k'] = k + 1
+    elif kind == 'radius':
+        params['k'] = features.shape[0]  # number of vertices
+        params['distance_upper_bound'] = radius
+    distances, neighbors = tree.query(features, **params)
+    if kind == 'knn':
+        return neighbors, distances
+    elif kind == 'radius':
+        dist = []
+        neigh = []
+        for distance, neighbor in zip(distances, neighbors):
+            mask = (distance != np.inf)
+            dist.append(distance[mask])
+            neigh.append(neighbor[mask])
+        return neigh, dist
+
+
+def _flann(features, metric, order, kind, k, radius, params):
+    if metric == 'max_dist':
+        raise ValueError('flann gives wrong results for metric="max_dist".')
+    try:
+        import cyflann as cfl
+    except Exception as e:
+        raise ImportError('Cannot import cyflann. Choose another nearest '
+                          'neighbors backend or try to install it with '
+                          'pip (or conda) install cyflann. '
+                          'Original exception: {}'.format(e))
+    cfl.set_distance_type(metric, order=order)
+    index = cfl.FLANNIndex()
+    index.build_index(features, **params)
+    # I tried changing the algorithm and testing performance on huge matrices,
+    # but the default parameters seems to work best.
+    if kind == 'knn':
+        neighbors, distances = index.nn_index(features, k+1)
+        if metric == 'euclidean':
+            np.sqrt(distances, out=distances)
+        elif metric == 'minkowski':
+            np.power(distances, 1/order, out=distances)
+    elif kind == 'radius':
+        distances = []
+        neighbors = []
+        if metric == 'euclidean':
+            radius = radius**2
+        elif metric == 'minkowski':
+            radius = radius**order
+        n_vertices, _ = features.shape
+        for vertex in range(n_vertices):
+            neighbor, distance = index.nn_radius(features[vertex, :], radius)
+            distances.append(distance)
+            neighbors.append(neighbor)
+        if metric == 'euclidean':
+            distances = list(map(np.sqrt, distances))
+        elif metric == 'minkowski':
+            distances = list(map(lambda d: np.power(d, 1/order), distances))
+    index.free_index()
+    return neighbors, distances
+
+
+def _nmslib(features, metric, order, kind, k, _, params):
+    if kind == 'radius':
+        raise ValueError('nmslib does not support kind="radius".')
+    if metric == 'minkowski':
+        raise ValueError('nmslib does not support metric="minkowski".')
+    try:
+        import nmslib as nms
+    except Exception as e:
+        raise ImportError('Cannot import nmslib. Choose another nearest '
+                          'neighbors backend or try to install it with '
+                          'pip (or conda) install nmslib. '
+                          'Original exception: {}'.format(e))
+    n_vertices, _ = features.shape
+    params_index = params.pop('index', None)
+    params_query = params.pop('query', None)
+    metric = 'l2' if metric == 'euclidean' else metric
+    metric = 'l1' if metric == 'manhattan' else metric
+    metric = 'linf' if metric == 'max_dist' else metric
+    index = nms.init(space=metric, **params)
+    index.addDataPointBatch(features)
+    index.createIndex(params_index)
+    if params_query is not None:
+        index.setQueryTimeParams(params_query)
+    results = index.knnQueryBatch(features, k=k+1)
+    neighbors, distances = zip(*results)
+    distances = np.concatenate(distances).reshape(n_vertices, k+1)
+    neighbors = np.concatenate(neighbors).reshape(n_vertices, k+1)
+    return neighbors, distances
+
+def nearest_neighbor(features, metric='euclidean', order=2, kind='knn', k=10, radius=None, backend='scipy-ckdtree', **kwargs):
+    '''Find nearest neighboors.
+    
+    Parameters
+    ----------
+    features : data numpy array 
+    metric : {'euclidean', 'manhattan', 'minkowski', 'max_dist'}, optional
+        Metric used to compute pairwise distances.
+
+        * ``'euclidean'`` defines pairwise distances as
+          :math:`d(v_i, v_j) = \| x_i - x_j \|_2`.
+        * ``'manhattan'`` defines pairwise distances as
+          :math:`d(v_i, v_j) = \| x_i - x_j \|_1`.
+        * ``'minkowski'`` generalizes the above and defines distances as
+          :math:`d(v_i, v_j) = \| x_i - x_j \|_p`
+          where :math:`p` is the ``order`` of the norm.
+        * ``'max_dist'`` defines pairwise distances as
+          :math:`d(v_i, v_j) = \| x_i - x_j \|_\infty = \max(x_i - x_j)`, where
+          the maximum is taken over the elements of the vector.
+
+        More metrics may be supported for some backends.
+        Please refer to the documentation of the chosen backend.
+    kind : 'knn' or 'radius' (default 'knn')
+    k : number of nearest neighboors if 'knn' is selected
+    radius : radius of the search if 'radius' is slected
+    
+    order : float, optional
+        The order of the Minkowski distance for ``metric='minkowski'``.
+    backend : string, optional
+        * ``'scipy-pdist'`` uses :func:`scipy.spatial.distance.pdist` to
+          compute pairwise distances. The method is brute force and computes
+          all distances. That is the slowest method.
+        * ``'scipy-kdtree'`` uses :class:`scipy.spatial.KDTree`. The method
+          builds a k-d tree to prune the number of pairwise distances it has to
+          compute. That is an efficient strategy for low-dimensional spaces.
+        * ``'scipy-ckdtree'`` uses :class:`scipy.spatial.cKDTree`. The same as
+          ``'scipy-kdtree'`` but with C bindings, which should be faster.
+          That is the default.
+        * ``'flann'`` uses the `Fast Library for Approximate Nearest Neighbors
+          (FLANN) <https://github.com/mariusmuja/flann>`_. That method is an
+          approximation.
+        * ``'nmslib'`` uses the `Non-Metric Space Library (NMSLIB)
+          <https://github.com/nmslib/nmslib>`_. That method is an
+          approximation. It should be the fastest in high-dimensional spaces.
+
+        You can look at this `benchmark
+        <https://github.com/erikbern/ann-benchmarks>`_ to get an idea of the
+        relative performance of those backends. It's nonetheless wise to run
+        some tests on your own data.   
+    '''
+    if kind=='knn':
+        radius = None
+    elif kind=='radius':
+        k = None
+    else:
+        raise ValueError('"kind" must be "knn" or "radius"')
+    
+    _orders = {
+        'euclidean': 2,
+        'manhattan': 1,
+        'max_dist': np.inf,
+        'minkowski': order,
+    }
+    order = _orders.pop(metric, None)  
+    try:
+        function = globals()['_' + backend.replace('-', '_')]
+    except KeyError:
+        raise ValueError('Invalid backend "{}".'.format(backend))
+    neighbors, distances = function(features, metric, order,
+                                    kind, k, radius, kwargs)
+    return neighbors, distances
+
+
+def sparse_distance_matrix(neighbors, distances, symmetrize=True, safe=False, kind = None, k=None):
+    '''Build a sparse distance matrix from nearest neighbors'''
+    n_edges = [len(n) - 1 for n in neighbors]  # remove distance to self
+    if safe and kind is None:
+        raise ValueError('Please specify "kind" to "knn" or "radius" to use the safe mode')
+    
+    n_vertices = len(n_edges)
+    if safe and kind == 'radius':
+        n_disconnected = np.sum(np.asarray(n_edges) == 0)
+        if n_disconnected > 0:
+            _logger.warning('{} points (out of {}) have no neighboors. '
+                            'Consider increasing the radius or setting '
+                            'kind=knn.'.format(n_disconnected, n_vertices))
+
+    value = np.empty(sum(n_edges), dtype=np.float)
+    row = np.empty_like(value, dtype=np.int)
+    col = np.empty_like(value, dtype=np.int)
+    start = 0
+    for vertex in range(n_vertices):
+        if safe and kind == 'knn':
+            assert n_edges[vertex] == k
+        end = start + n_edges[vertex]
+        value[start:end] = distances[vertex][1:]
+        row[start:end] = np.full(n_edges[vertex], vertex)
+        col[start:end] = neighbors[vertex][1:]
+        start = end
+    W = sparse.csr_matrix((value, (row, col)), (n_vertices, n_vertices))
+    if symmetrize:
+        # Enforce symmetry. May have been broken by k-NN. Checking symmetry
+        # with np.abs(W - W.T).sum() is as costly as the symmetrization itself.
+        W = utils.symmetrize(W, method='fill')
+    return W
\ No newline at end of file
diff --git a/pygsp/filters/filter.py b/pygsp/filters/filter.py
index fca48b37..7e36f11c 100644
--- a/pygsp/filters/filter.py
+++ b/pygsp/filters/filter.py
@@ -255,7 +255,7 @@ def filter(self, s, method='chebyshev', order=30):
         >>> _ = G.plot(s1, ax=axes[0])
         >>> _ = G.plot(s2, ax=axes[1])
         >>> print('{:.5f}'.format(np.linalg.norm(s1 - s2)))
-        0.26808
+        0.26995
 
         Perfect reconstruction with Itersine, a tight frame:
 
@@ -448,7 +448,7 @@ def estimate_frame_bounds(self, x=None):
         A=1.708, B=2.359
         >>> A, B = g.estimate_frame_bounds(G.e)
         >>> print('A={:.3f}, B={:.3f}'.format(A, B))
-        A=1.723, B=2.359
+        A=1.839, B=2.359
 
         The frame bounds can be seen in the plot of the filter bank as the
         minimum and maximum of their squared sum (the black curve):
diff --git a/pygsp/graphs/fourier.py b/pygsp/graphs/fourier.py
index 3c215273..d699ceaa 100644
--- a/pygsp/graphs/fourier.py
+++ b/pygsp/graphs/fourier.py
@@ -85,7 +85,7 @@ def coherence(self):
         >>> graph.compute_fourier_basis()
         >>> minimum = 1 / np.sqrt(graph.n_vertices)
         >>> print('{:.2f} in [{:.2f}, 1]'.format(graph.coherence, minimum))
-        0.75 in [0.12, 1]
+        0.87 in [0.12, 1]
         >>>
         >>> # Plot the most localized eigenvector.
         >>> import matplotlib.pyplot as plt
diff --git a/pygsp/graphs/nngraphs/bunny.py b/pygsp/graphs/nngraphs/bunny.py
index d9dac407..cc269c10 100644
--- a/pygsp/graphs/nngraphs/bunny.py
+++ b/pygsp/graphs/nngraphs/bunny.py
@@ -34,7 +34,5 @@ def __init__(self, **kwargs):
             'distance': 8,
         }
 
-        super(Bunny, self).__init__(Xin=data['bunny'],
-                                    epsilon=0.02, NNtype='radius',
-                                    center=False, rescale=False,
+        super(Bunny, self).__init__(data['bunny'], kind='radius', radius=0.02,
                                     plotting=plotting, **kwargs)
diff --git a/pygsp/graphs/nngraphs/cube.py b/pygsp/graphs/nngraphs/cube.py
index 820e401c..b9cb1d60 100644
--- a/pygsp/graphs/nngraphs/cube.py
+++ b/pygsp/graphs/nngraphs/cube.py
@@ -10,12 +10,12 @@ class Cube(NNGraph):
 
     Parameters
     ----------
-    radius : float
-        Edge lenght (default = 1)
     nb_pts : int
         Number of vertices (default = 300)
     nb_dim : int
         Dimension (default = 3)
+    length : float
+        Edge length (default = 1)
     sampling : string
         Variance of the distance kernel (default = 'random')
         (Can now only be 'random')
@@ -35,20 +35,23 @@ class Cube(NNGraph):
     """
 
     def __init__(self,
-                 radius=1,
                  nb_pts=300,
                  nb_dim=3,
+                 length=1,
                  sampling='random',
                  seed=None,
                  **kwargs):
 
-        self.radius = radius
         self.nb_pts = nb_pts
         self.nb_dim = nb_dim
+        self.length = length
         self.sampling = sampling
         self.seed = seed
         rs = np.random.RandomState(seed)
 
+        if length != 1:
+            raise NotImplementedError('Only length=1 is implemented.')
+
         if self.nb_dim > 3:
             raise NotImplementedError("Dimension > 3 not supported yet!")
 
@@ -89,12 +92,10 @@ def __init__(self,
             'distance': 9,
         }
 
-        super(Cube, self).__init__(Xin=pts, k=10,
-                                   center=False, rescale=False,
-                                   plotting=plotting, **kwargs)
+        super(Cube, self).__init__(pts, k=10, plotting=plotting, **kwargs)
 
     def _get_extra_repr(self):
-        attrs = {'radius': '{:.2f}'.format(self.radius),
+        attrs = {'length': '{:.2e}'.format(self.length),
                  'nb_pts': self.nb_pts,
                  'nb_dim': self.nb_dim,
                  'sampling': self.sampling,
diff --git a/pygsp/graphs/nngraphs/imgpatches.py b/pygsp/graphs/nngraphs/imgpatches.py
index 73cb6abd..df681a3b 100644
--- a/pygsp/graphs/nngraphs/imgpatches.py
+++ b/pygsp/graphs/nngraphs/imgpatches.py
@@ -10,7 +10,8 @@ class ImgPatches(NNGraph):
 
     Extract a feature vector in the form of a patch for every pixel of an
     image, then construct a nearest-neighbor graph between these feature
-    vectors. The feature matrix, i.e. the patches, can be found in :attr:`Xin`.
+    vectors. The feature matrix, i.e., the patches, can be found in
+    :attr:`features`.
 
     Parameters
     ----------
@@ -39,9 +40,10 @@ class ImgPatches(NNGraph):
     >>> from skimage import data, img_as_float
     >>> img = img_as_float(data.camera()[::64, ::64])
     >>> G = graphs.ImgPatches(img, patch_shape=(3, 3))
-    >>> print('{} nodes ({} x {} pixels)'.format(G.Xin.shape[0], *img.shape))
+    >>> N, d = G.patches.shape
+    >>> print('{} nodes ({} x {} pixels)'.format(N, *img.shape))
     64 nodes (8 x 8 pixels)
-    >>> print('{} features per node'.format(G.Xin.shape[1]))
+    >>> print('{} features per node'.format(d))
     9 features per node
     >>> G.set_coordinates(kind='spring', seed=42)
     >>> fig, axes = plt.subplots(1, 2)
@@ -87,16 +89,16 @@ def __init__(self, img, patch_shape=(3, 3), **kwargs):
         # Alternative: sklearn.feature_extraction.image.extract_patches_2d.
         #              sklearn has much less dependencies than skimage.
         try:
-            import skimage
+            from skimage.util import view_as_windows
         except Exception as e:
             raise ImportError('Cannot import skimage, which is needed to '
                               'extract patches. Try to install it with '
                               'pip (or conda) install scikit-image. '
                               'Original exception: {}'.format(e))
-        patches = skimage.util.view_as_windows(img, window_shape=window_shape)
-        patches = patches.reshape((h * w, r * c * d))
+        self.patches = view_as_windows(img, window_shape)
+        self.patches = self.patches.reshape((h * w, r * c * d))
 
-        super(ImgPatches, self).__init__(patches, **kwargs)
+        super(ImgPatches, self).__init__(self.patches, **kwargs)
 
     def _get_extra_repr(self):
         attrs = dict(patch_shape=self.patch_shape)
diff --git a/pygsp/graphs/nngraphs/nngraph.py b/pygsp/graphs/nngraphs/nngraph.py
index 252a5416..3f56734e 100644
--- a/pygsp/graphs/nngraphs/nngraph.py
+++ b/pygsp/graphs/nngraphs/nngraph.py
@@ -1,199 +1,401 @@
 # -*- coding: utf-8 -*-
 
-import traceback
+from __future__ import division
+
+from functools import partial
 
 import numpy as np
 from scipy import sparse, spatial
 
 from pygsp import utils
 from pygsp.graphs import Graph  # prevent circular import in Python < 3.5
+from pygsp._nearest_neighbor import nearest_neighbor, sparse_distance_matrix
+
 
 _logger = utils.build_logger(__name__)
 
 
-def _import_pfl():
-    try:
-        import pyflann as pfl
-    except Exception as e:
-        raise ImportError('Cannot import pyflann. Choose another nearest '
-                          'neighbors method or try to install it with '
-                          'pip (or conda) install pyflann (or pyflann3). '
-                          'Original exception: {}'.format(e))
-    return pfl
+class NNGraph(Graph):
+    r"""Nearest-neighbor graph.
+
+    The nearest-neighbor graph is built from a set of features, where the edge
+    weight between vertices :math:`v_i` and :math:`v_j` is given by
 
+    .. math:: A(i,j) = k \left( \frac{d(v_i, v_j)}{\sigma} \right),
 
-class NNGraph(Graph):
-    r"""Nearest-neighbor graph from given point cloud.
+    where :math:`d(v_i, v_j)` is a distance measure between some representation
+    (the features) of :math:`v_i` and :math:`v_j`, :math:`k` is a kernel
+    function that transforms a distance in :math:`[0, \infty]` to a similarity
+    measure generally in :math:`[0, 1]`, and :math:`\sigma` is the kernel width.
+
+    For example, the features might be the 3D coordinates of points in a point
+    cloud. Then, if ``metric='euclidean'`` and ``kernel='gaussian'`` (the
+    defaults), :math:`A(i,j) = \exp(-\log(2) \| x_i - x_j \|_2^2 / \sigma^2)`,
+    where :math:`x_i` is the 3D position of vertex :math:`v_i`.
+
+    The similarity matrix :math:`A` is sparsified by either keeping the ``k``
+    closest vertices for each vertex (if ``type='knn'``), or by setting to zero
+    the similarity when the distance is greater than ``radius`` (if ``type='radius'``).
 
     Parameters
     ----------
-    Xin : ndarray
-        Input points, Should be an `N`-by-`d` matrix, where `N` is the number
-        of nodes in the graph and `d` is the dimension of the feature space.
-    NNtype : string, optional
-        Type of nearest neighbor graph to create. The options are 'knn' for
-        k-Nearest Neighbors or 'radius' for epsilon-Nearest Neighbors (default
-        is 'knn').
-    use_flann : bool, optional
-        Use Fast Library for Approximate Nearest Neighbors (FLANN) or not.
-        (default is False)
-    center : bool, optional
-        Center the data so that it has zero mean (default is True)
-    rescale : bool, optional
-        Rescale the data so that it lies in a l2-sphere (default is True)
-    k : int, optional
-        Number of neighbors for knn (default is 10)
-    sigma : float, optional
-        Width of the similarity kernel.
-        By default, it is set to the average of the nearest neighbor distance.
-    epsilon : float, optional
-        Radius for the epsilon-neighborhood search (default is 0.01)
-    plotting : dict, optional
-        Dictionary of plotting parameters. See :obj:`pygsp.plotting`.
-        (default is {})
-    symmetrize_type : string, optional
-        Type of symmetrization to use for the adjacency matrix. See
-        :func:`pygsp.utils.symmetrization` for the options.
-        (default is 'average')
-    dist_type : string, optional
-        Type of distance to compute. See
-        :func:`pyflann.index.set_distance_type` for possible options.
-        (default is 'euclidean')
+    features : ndarray
+        An `N`-by-`d` matrix, where `N` is the number of nodes in the graph and
+        `d` is the number of features.
+    standardize : bool, optional
+        Whether to rescale the features so that each feature has a mean of 0
+        and standard deviation of 1 (unit variance).
+    metric : {'euclidean', 'manhattan', 'minkowski', 'max_dist'}, optional
+        Metric used to compute pairwise distances.
+
+        * ``'euclidean'`` defines pairwise distances as
+          :math:`d(v_i, v_j) = \| x_i - x_j \|_2`.
+        * ``'manhattan'`` defines pairwise distances as
+          :math:`d(v_i, v_j) = \| x_i - x_j \|_1`.
+        * ``'minkowski'`` generalizes the above and defines distances as
+          :math:`d(v_i, v_j) = \| x_i - x_j \|_p`
+          where :math:`p` is the ``order`` of the norm.
+        * ``'max_dist'`` defines pairwise distances as
+          :math:`d(v_i, v_j) = \| x_i - x_j \|_\infty = \max(x_i - x_j)`, where
+          the maximum is taken over the elements of the vector.
+
+        More metrics may be supported for some backends.
+        Please refer to the documentation of the chosen backend.
     order : float, optional
-        Only used if dist_type is 'minkowski'; represents the order of the
-        Minkowski distance. (default is 0)
+        The order of the Minkowski distance for ``metric='minkowski'``.
+    kind : {'knn', 'radius'}, optional
+        Kind of nearest neighbor graph to create. Either ``'knn'`` for
+        k-nearest neighbors or ``'radius'`` for epsilon-nearest neighbors.
+    k : int, optional
+        Number of neighbors considered when building a k-NN graph with
+        ``type='knn'``.
+    radius : float or {'estimate', 'estimate-knn'}, optional
+        Radius of the ball when building a radius graph with ``type='radius'``.
+        It is hard to set an optimal radius. If too small, some vertices won't
+        be connected to any other vertex. If too high, vertices will be
+        connected to many other vertices and the graph won't be sparse (high
+        average degree).  If no good radius is known a priori, we can estimate
+        one. ``'estimate'`` sets the radius as the expected average distance
+        between vertices for a uniform sampling of the ambient space.
+        ``'estimate-knn'`` first builds a knn graph and sets the radius to the
+        average distance. ``'estimate-knn'`` usually gives a better estimation
+        but is more costly. ``'estimate'`` can be better in low dimension.
+    kernel : string or function
+        The function :math:`k` that transforms a distance to a similarity.
+        The following kernels are pre-defined.
+
+        * ``'gaussian'`` defines the Gaussian, also known as the radial basis
+          function (RBF), kernel :math:`k(d) = \exp(-\log(2) d^2)`.
+        * ``'exponential'`` defines the kernel :math:`k(d) = \exp(-\log(2) d)`.
+        * ``'rectangular'`` returns 1 if :math:`d < 1` and 0 otherwise.
+        * ``'triangular'`` defines the kernel :math:`k(d) = \max(1 - d/2, 0)`.
+        * Other kernels are ``'tricube'``, ``'triweight'``, ``'quartic'``,
+          ``'epanechnikov'``, ``'logistic'``, and ``'sigmoid'``.
+          See `Wikipedia <https://en.wikipedia.org/wiki/Kernel_(statistics)>`_.
+
+        Another option is to pass a function that takes a vector of pairwise
+        distances and returns the similarities. All the predefined kernels
+        return a similarity of 0.5 when the distance is one.
+        An example of custom kernel is ``kernel=lambda d: d.min() / d``.
+    kernel_width : float, optional
+        Control the width, also known as the bandwidth, :math:`\sigma` of the
+        kernel. It scales the distances as ``distances / kernel_width`` before
+        calling the kernel function.
+        By default, it is set to the average of all computed distances for
+        ``kind='knn'`` and to half the radius for ``kind='radius'``.
+    backend : string, optional
+        * ``'scipy-pdist'`` uses :func:`scipy.spatial.distance.pdist` to
+          compute pairwise distances. The method is brute force and computes
+          all distances. That is the slowest method.
+        * ``'scipy-kdtree'`` uses :class:`scipy.spatial.KDTree`. The method
+          builds a k-d tree to prune the number of pairwise distances it has to
+          compute. That is an efficient strategy for low-dimensional spaces.
+        * ``'scipy-ckdtree'`` uses :class:`scipy.spatial.cKDTree`. The same as
+          ``'scipy-kdtree'`` but with C bindings, which should be faster.
+          That is the default.
+        * ``'flann'`` uses the `Fast Library for Approximate Nearest Neighbors
+          (FLANN) <https://github.com/mariusmuja/flann>`_. That method is an
+          approximation.
+        * ``'nmslib'`` uses the `Non-Metric Space Library (NMSLIB)
+          <https://github.com/nmslib/nmslib>`_. That method is an
+          approximation. It should be the fastest in high-dimensional spaces.
+
+        You can look at this `benchmark
+        <https://github.com/erikbern/ann-benchmarks>`_ to get an idea of the
+        relative performance of those backends. It's nonetheless wise to run
+        some tests on your own data.
+    kwargs : dict
+        Parameters to be passed to the :class:`Graph` constructor or the
+        backend library.
 
     Examples
     --------
+
+    Construction of a graph from a set of features.
+
     >>> import matplotlib.pyplot as plt
-    >>> X = np.random.RandomState(42).uniform(size=(30, 2))
-    >>> G = graphs.NNGraph(X)
+    >>> rs = np.random.RandomState(42)
+    >>> features = rs.uniform(size=(30, 2))
+    >>> G = graphs.NNGraph(features)
     >>> fig, axes = plt.subplots(1, 2)
     >>> _ = axes[0].spy(G.W, markersize=5)
     >>> _ = G.plot(ax=axes[1])
 
-    """
+    Radius versus knn graph.
+
+    >>> features = rs.uniform(size=(100, 3))
+    >>> fig, ax = plt.subplots()
+    >>> G = graphs.NNGraph(features, kind='radius', radius=0.2964)
+    >>> label = 'radius graph ({} edges)'.format(G.n_edges)
+    >>> _ = ax.hist(G.W.data, bins=20, label=label, alpha=0.5)
+    >>> G = graphs.NNGraph(features, kind='knn', k=6)
+    >>> label = 'knn graph ({} edges)'.format(G.n_edges)
+    >>> _ = ax.hist(G.W.data, bins=20, label=label, alpha=0.5)
+    >>> _ = ax.legend()
+    >>> _ = ax.set_title('edge weights')
+
+    Control of the sparsity of knn and radius graphs.
+
+    >>> features = rs.uniform(size=(100, 3))
+    >>> n_edges = dict(knn=[], radius=[])
+    >>> n_neighbors = np.arange(1, 100, 5)
+    >>> radiuses = np.arange(0.05, 1.5, 0.05)
+    >>> for k in n_neighbors:
+    ...     G = graphs.NNGraph(features, kind='knn', k=k)
+    ...     n_edges['knn'].append(G.n_edges)
+    >>> for radius in radiuses:
+    ...     G = graphs.NNGraph(features, kind='radius', radius=radius)
+    ...     n_edges['radius'].append(G.n_edges)
+    >>> fig, axes = plt.subplots(1, 2, sharey=True)
+    >>> _ = axes[0].plot(n_neighbors, n_edges['knn'])
+    >>> _ = axes[1].plot(radiuses, n_edges['radius'])
+    >>> _ = axes[0].set_ylabel('number of edges')
+    >>> _ = axes[0].set_xlabel('number of neighbors (knn graph)')
+    >>> _ = axes[1].set_xlabel('radius (radius graph)')
+    >>> _ = fig.suptitle('Sparsity')
+
+    Choice of metric and the curse of dimensionality.
 
-    def __init__(self, Xin, NNtype='knn', use_flann=False, center=True,
-                 rescale=True, k=10, sigma=None, epsilon=0.01,
-                 plotting={}, symmetrize_type='average', dist_type='euclidean',
-                 order=0, **kwargs):
+    >>> fig, axes = plt.subplots(1, 2)
+    >>> for dim, ax in zip([3, 30], axes):
+    ...     features = rs.uniform(size=(100, dim))
+    ...     for metric in ['euclidean', 'manhattan', 'max_dist', 'cosine']:
+    ...         G = graphs.NNGraph(features, metric=metric,
+    ...                            backend='scipy-pdist')
+    ...         _ = ax.hist(G.W.data, bins=20, label=metric, alpha=0.5)
+    ...     _ = ax.legend()
+    ...     _ = ax.set_title('edge weights, {} dimensions'.format(dim))
 
-        self.Xin = Xin
-        self.NNtype = NNtype
-        self.use_flann = use_flann
-        self.center = center
-        self.rescale = rescale
-        self.k = k
-        self.sigma = sigma
-        self.epsilon = epsilon
-        self.symmetrize_type = symmetrize_type
-        self.dist_type = dist_type
-        self.order = order
+    Choice of kernel.
 
-        N, d = np.shape(self.Xin)
-        Xout = self.Xin
-
-        if k >= N:
-            raise ValueError('The number of neighbors (k={}) must be smaller '
-                             'than the number of nodes ({}).'.format(k, N))
-
-        if self.center:
-            Xout = self.Xin - np.kron(np.ones((N, 1)),
-                                      np.mean(self.Xin, axis=0))
-
-        if self.rescale:
-            bounding_radius = 0.5 * np.linalg.norm(np.amax(Xout, axis=0) -
-                                                   np.amin(Xout, axis=0), 2)
-            scale = np.power(N, 1. / float(min(d, 3))) / 10.
-            Xout *= scale / bounding_radius
-
-        # Translate distance type string to corresponding Minkowski order.
-        dist_translation = {"euclidean": 2,
-                            "manhattan": 1,
-                            "max_dist": np.inf,
-                            "minkowski": order
-                            }
-
-        if self.NNtype == 'knn':
-            spi = np.zeros((N * k))
-            spj = np.zeros((N * k))
-            spv = np.zeros((N * k))
-
-            if self.use_flann:
-                pfl = _import_pfl()
-                pfl.set_distance_type(dist_type, order=order)
-                flann = pfl.FLANN()
-
-                # Default FLANN parameters (I tried changing the algorithm and
-                # testing performance on huge matrices, but the default one
-                # seems to work best).
-                NN, D = flann.nn(Xout, Xout, num_neighbors=(k + 1),
-                                 algorithm='kdtree')
-
-            else:
-                kdt = spatial.KDTree(Xout)
-                D, NN = kdt.query(Xout, k=(k + 1),
-                                  p=dist_translation[dist_type])
-
-            if self.sigma is None:
-                self.sigma = np.mean(D[:, 1:])  # Discard distance to self.
-
-            for i in range(N):
-                spi[i * k:(i + 1) * k] = np.kron(np.ones((k)), i)
-                spj[i * k:(i + 1) * k] = NN[i, 1:]
-                spv[i * k:(i + 1) * k] = np.exp(-np.power(D[i, 1:], 2) /
-                                                float(self.sigma))
-
-        elif self.NNtype == 'radius':
-
-            kdt = spatial.KDTree(Xout)
-            D, NN = kdt.query(Xout, k=None, distance_upper_bound=epsilon,
-                              p=dist_translation[dist_type])
-            if self.sigma is None:
-                # Discard distance to self.
-                self.sigma = np.mean([np.mean(d[1:]) for d in D])
-            count = 0
-            for i in range(N):
-                count = count + len(NN[i])
-
-            spi = np.zeros((count))
-            spj = np.zeros((count))
-            spv = np.zeros((count))
-
-            start = 0
-            for i in range(N):
-                leng = len(NN[i]) - 1
-                spi[start:start + leng] = np.kron(np.ones((leng)), i)
-                spj[start:start + leng] = NN[i][1:]
-                spv[start:start + leng] = np.exp(-np.power(D[i][1:], 2) /
-                                                 float(self.sigma))
-                start = start + leng
+    >>> fig, axes = plt.subplots(1, 2)
+    >>> width = 0.3
+    >>> distances = np.linspace(0, 1, 200)
+    >>> for name, kernel in graphs.NNGraph._kernels.items():
+    ...     _ = axes[0].plot(distances, kernel(distances / width), label=name)
+    >>> _ = axes[0].set_xlabel('distance [0, inf]')
+    >>> _ = axes[0].set_ylabel('similarity [0, 1]')
+    >>> _ = axes[0].legend(loc='upper right')
+    >>> features = rs.uniform(size=(100, 3))
+    >>> for kernel in ['gaussian', 'triangular', 'tricube', 'exponential']:
+    ...     G = graphs.NNGraph(features, kernel=kernel)
+    ...     _ = axes[1].hist(G.W.data, bins=20, label=kernel, alpha=0.5)
+    >>> _ = axes[1].legend()
+    >>> _ = axes[1].set_title('edge weights')
+
+    Choice of kernel width.
+
+    >>> fig, axes = plt.subplots()
+    >>> for width in [.2, .3, .4, .6, .8, None]:
+    ...     G = graphs.NNGraph(features, kernel_width=width)
+    ...     label = 'width = {:.2f}'.format(G.kernel_width)
+    ...     _ = axes.hist(G.W.data, bins=20, label=label, alpha=0.5)
+    >>> _ = axes.legend(loc='upper left')
+    >>> _ = axes.set_title('edge weights')
+
+    Choice of backend. Compare on your data!
+
+    >>> import time
+    >>> sizes = [300, 1000, 3000]
+    >>> dims = [3, 100]
+    >>> backends = ['scipy-pdist', 'scipy-kdtree', 'scipy-ckdtree', 'flann',
+    ...             'nmslib']
+    >>> times = np.full((len(sizes), len(dims), len(backends)), np.nan)
+    >>> for i, size in enumerate(sizes):
+    ...     for j, dim in enumerate(dims):
+    ...         for k, backend in enumerate(backends):
+    ...             if (size * dim) > 1e4 and backend == 'scipy-kdtree':
+    ...                 continue  # too slow
+    ...             features = rs.uniform(size=(size, dim))
+    ...             start = time.time()
+    ...             _ = graphs.NNGraph(features, backend=backend)
+    ...             times[i][j][k] = time.time() - start
+    >>> fig, axes = plt.subplots(1, 2, sharey=True)
+    >>> for j, (dim, ax) in enumerate(zip(dims, axes)):
+    ...     for k, backend in enumerate(backends):
+    ...         _ = ax.loglog(sizes, times[:, j, k], '.-', label=backend)
+    ...         _ = ax.set_title('{} dimensions'.format(dim))
+    ...         _ = ax.set_xlabel('number of vertices')
+    >>> _ = axes[0].set_ylabel('execution time [s]')
+    >>> _ = axes[1].legend(loc='upper left')
 
+    """
+
+    def __init__(self, features, standardize=False,
+                 metric='euclidean', order=3,
+                 kind='knn', k=10, radius='estimate-knn',
+                 kernel='gaussian', kernel_width=None,
+                 backend='scipy-ckdtree',
+                 **kwargs):
+
+        # features is stored in coords, potentially standardized
+        self.standardize = standardize
+        self.metric = metric
+        self.kind = kind
+        self.kernel = kernel
+        self.k = k
+        self.backend = backend
+
+        features = np.asanyarray(features)
+        if features.ndim != 2:
+            raise ValueError('features should be #vertices x dimensionality')
+        n_vertices, dimensionality = features.shape
+
+        params_graph = dict()
+        for key in ['lap_type', 'plotting']:
+            try:
+                params_graph[key] = kwargs.pop(key)
+            except KeyError:
+                pass
+
+        if standardize:
+            # Don't alter the original data (users would be surprised).
+            features = features - np.mean(features, axis=0)
+            features /= np.std(features, axis=0)
+
+        # Order consistent with metric (used by kdtree and ckdtree).
+        _orders = {
+            'euclidean': 2,
+            'manhattan': 1,
+            'max_dist': np.inf,
+            'minkowski': order,
+        }
+        order = _orders.pop(metric, None)
+
+        if kind == 'knn':
+            if not 1 <= k < n_vertices:
+                raise ValueError('The number of neighbors (k={}) must be '
+                                 'greater than 0 and smaller than the number '
+                                 'of vertices ({}).'.format(k, n_vertices))
+            radius = None
+        elif kind == 'radius':
+            if radius == 'estimate':
+                maximums = np.amax(features, axis=0)
+                minimums = np.amin(features, axis=0)
+                distance_max = np.linalg.norm(maximums - minimums, order)
+                radius = distance_max / np.power(n_vertices, 1/dimensionality)
+            elif radius == 'estimate-knn':
+                graph = NNGraph(features, standardize=standardize,
+                                metric=metric, order=order, kind='knn', k=k,
+                                kernel_width=None, backend=backend, **kwargs)
+                radius = graph.kernel_width
+            elif radius <= 0:
+                raise ValueError('The radius must be greater than 0.')
+            self.k = None
         else:
-            raise ValueError('Unknown NNtype {}'.format(self.NNtype))
+            raise ValueError('Invalid kind "{}".'.format(kind))
+
+        neighbors, distances = nearest_neighbor(features, metric=metric, order=order,
+                                        kind=kind, k=k, radius=radius, backend=backend, **kwargs)
+        W = sparse_distance_matrix(neighbors, distances, symmetrize=True, safe=True, kind = kind, k=k)
 
-        W = sparse.csc_matrix((spv, (spi, spj)), shape=(N, N))
+        if kernel_width is None:
+            if kind == 'knn':
+                kernel_width = np.mean(W.data) if W.nnz > 0 else np.nan
+            elif kind == 'radius':
+                kernel_width = radius / 2
 
-        # Sanity check
-        if np.shape(W)[0] != np.shape(W)[1]:
-            raise ValueError('Weight matrix W is not square')
+        if not callable(kernel):
+            try:
+                kernel = self._kernels[kernel]
+            except KeyError:
+                raise ValueError('Unknown kernel {}.'.format(kernel))
 
-        # Enforce symmetry. Note that checking symmetry with
-        # np.abs(W - W.T).sum() is as costly as the symmetrization itself.
-        W = utils.symmetrize(W, method=symmetrize_type)
+        assert np.all(W.data >= 0), 'Distance must be in [0, inf].'
+        W.data = kernel(W.data / kernel_width)
+        if not np.all((W.data >= 0) & (W.data <= 1)):
+            _logger.warning('Kernel returned similarity not in [0, 1].')
+
+        self.order = order
+        self.radius = radius
+        self.kernel_width = kernel_width
 
-        super(NNGraph, self).__init__(W, plotting=plotting,
-                                      coords=Xout, **kwargs)
+        super(NNGraph, self).__init__(W, coords=features, **params_graph)
 
     def _get_extra_repr(self):
-        return {'NNtype': self.NNtype,
-                'use_flann': self.use_flann,
-                'center': self.center,
-                'rescale': self.rescale,
-                'k': self.k,
-                'sigma': '{:.2f}'.format(self.sigma),
-                'epsilon': '{:.2f}'.format(self.epsilon),
-                'symmetrize_type': self.symmetrize_type,
-                'dist_type': self.dist_type,
-                'order': self.order}
+        attrs = {
+            'standardize': self.standardize,
+            'metric': self.metric,
+            'order': self.order,
+            'kind': self.kind,
+        }
+        if self.k is not None:
+            attrs['k'] = self.k
+        if self.radius is not None:
+            attrs['radius'] = '{:.2e}'.format(self.radius)
+        attrs.update({
+            'kernel': '{}'.format(self.kernel),
+            'kernel_width': '{:.2e}'.format(self.kernel_width),
+            'backend': self.backend,
+        })
+        return attrs
+
+    @staticmethod
+    def _kernel_rectangular(distance):
+        return (distance < 1).astype(np.float)
+
+    @staticmethod
+    def _kernel_triangular(distance, value_at_one=0.5):
+        distance = value_at_one * distance
+        return np.maximum(1 - distance, 0)
+
+    @staticmethod
+    def _kernel_exponential(distance, power=1, value_at_one=0.5):
+        cst = np.log(value_at_one)
+        return np.exp(cst * distance**power)
+
+    @staticmethod
+    def _kernel_powers(distance, pow1, pow2, value_at_one=0.5):
+        cst = (1 - value_at_one**(1/pow2))**(1/pow1)
+        distance = np.clip(cst * distance, 0, 1)
+        return (1 - distance**pow1)**pow2
+
+    @staticmethod
+    def _kernel_logistic(distance, value_at_one=0.5):
+        cst = 4 / value_at_one - 2
+        cst = np.log(0.5 * (cst + np.sqrt(cst**2 - 4)))
+        distance = cst * distance
+        return 4 / (np.exp(distance) + 2 + np.exp(-distance))
+
+    @staticmethod
+    def _kernel_sigmoid(distance, value_at_one=0.5):
+        cst = 2 / value_at_one
+        cst = np.log(0.5 * (cst + np.sqrt(cst**2 - 4)))
+        distance = cst * distance
+        return 2 / (np.exp(distance) + np.exp(-distance))
+
+    _kernels = {
+        'rectangular': _kernel_rectangular.__func__,
+        'triangular': _kernel_triangular.__func__,
+
+        'exponential': _kernel_exponential.__func__,
+        'gaussian': partial(_kernel_exponential.__func__, power=2),
+
+        'tricube': partial(_kernel_powers.__func__, pow1=3, pow2=3),
+        'triweight': partial(_kernel_powers.__func__, pow1=2, pow2=3),
+        'quartic': partial(_kernel_powers.__func__, pow1=2, pow2=2),
+        'epanechnikov': partial(_kernel_powers.__func__, pow1=2, pow2=1),
+
+        'logistic': _kernel_logistic.__func__,
+        'sigmoid': _kernel_sigmoid.__func__,
+    }
diff --git a/pygsp/graphs/nngraphs/sensor.py b/pygsp/graphs/nngraphs/sensor.py
index ac623a50..09764b02 100644
--- a/pygsp/graphs/nngraphs/sensor.py
+++ b/pygsp/graphs/nngraphs/sensor.py
@@ -75,9 +75,7 @@ def __init__(self, N=64, k=6, distributed=False, seed=None, **kwargs):
 
             coords = rs.uniform(0, 1, (N, 2))
 
-        super(Sensor, self).__init__(Xin=coords, k=k,
-                                     rescale=False, center=False,
-                                     plotting=plotting, **kwargs)
+        super(Sensor, self).__init__(coords, k=k, plotting=plotting, **kwargs)
 
     def _get_extra_repr(self):
         return {'k': self.k,
diff --git a/pygsp/graphs/nngraphs/sphere.py b/pygsp/graphs/nngraphs/sphere.py
index e375b5b4..1262b922 100644
--- a/pygsp/graphs/nngraphs/sphere.py
+++ b/pygsp/graphs/nngraphs/sphere.py
@@ -10,12 +10,12 @@ class Sphere(NNGraph):
 
     Parameters
     ----------
-    radius : float
-        Radius of the sphere (default = 1)
     nb_pts : int
         Number of vertices (default = 300)
     nb_dim : int
         Dimension (default = 3)
+    diameter : float
+        Radius of the sphere (default = 2)
     sampling : string
         Variance of the distance kernel (default = 'random')
         (Can now only be 'random')
@@ -35,14 +35,14 @@ class Sphere(NNGraph):
     """
 
     def __init__(self,
-                 radius=1,
                  nb_pts=300,
                  nb_dim=3,
+                 diameter=2,
                  sampling='random',
                  seed=None,
                  **kwargs):
 
-        self.radius = radius
+        self.diameter = diameter
         self.nb_pts = nb_pts
         self.nb_dim = nb_dim
         self.sampling = sampling
@@ -55,6 +55,7 @@ def __init__(self,
 
             for i in range(self.nb_pts):
                 pts[i] /= np.linalg.norm(pts[i])
+                pts[i] *= (diameter / 2)
 
         else:
 
@@ -64,12 +65,10 @@ def __init__(self,
             'vertex_size': 80,
         }
 
-        super(Sphere, self).__init__(Xin=pts, k=10,
-                                     center=False, rescale=False,
-                                     plotting=plotting, **kwargs)
+        super(Sphere, self).__init__(pts, k=10, plotting=plotting, **kwargs)
 
     def _get_extra_repr(self):
-        attrs = {'radius': '{:.2f}'.format(self.radius),
+        attrs = {'diameter': '{:.2e}'.format(self.diameter),
                  'nb_pts': self.nb_pts,
                  'nb_dim': self.nb_dim,
                  'sampling': self.sampling,
diff --git a/pygsp/graphs/nngraphs/twomoons.py b/pygsp/graphs/nngraphs/twomoons.py
index af87d0d8..1d7cdde1 100644
--- a/pygsp/graphs/nngraphs/twomoons.py
+++ b/pygsp/graphs/nngraphs/twomoons.py
@@ -16,7 +16,7 @@ class TwoMoons(NNGraph):
         two_moons graph or a synthesized one (default is 'standard').
         'standard' : Create a two_moons graph from a based graph.
         'synthesized' : Create a synthesized two_moon
-    sigmag : float
+    kernel_width : float
         Variance of the distance kernel (default = 0.05)
     dim : int
         The dimensionality of the points (default = 2).
@@ -24,7 +24,7 @@ class TwoMoons(NNGraph):
     N : int
         Number of vertices (default = 2000)
         Only valid for moontype == 'synthesized'.
-    sigmad : float
+    variance : float
         Variance of the data (do not set it too high or you won't see anything)
         (default = 0.05)
         Only valid for moontype == 'synthesized'.
@@ -44,11 +44,11 @@ class TwoMoons(NNGraph):
 
     """
 
-    def _create_arc_moon(self, N, sigmad, distance, number, seed):
+    def _create_arc_moon(self, N, variance, distance, number, seed):
         rs = np.random.RandomState(seed)
         phi = rs.rand(N, 1) * np.pi
         r = 1
-        rb = sigmad * rs.normal(size=(N, 1))
+        rb = variance * rs.normal(size=(N, 1))
         ab = rs.rand(N, 1) * 2 * np.pi
         b = rb * np.exp(1j * ab)
         bx = np.real(b)
@@ -63,29 +63,29 @@ def _create_arc_moon(self, N, sigmad, distance, number, seed):
 
         return np.concatenate((moonx, moony), axis=1)
 
-    def __init__(self, moontype='standard', dim=2, sigmag=0.05,
-                 N=400, sigmad=0.07, distance=0.5, seed=None, **kwargs):
+    def __init__(self, moontype='standard', dim=2, kernel_width=0.05,
+                 N=400, variance=0.07, distance=0.5, seed=None, **kwargs):
 
         self.moontype = moontype
         self.dim = dim
-        self.sigmag = sigmag
-        self.sigmad = sigmad
+        self.kernel_width = kernel_width
+        self.variance = variance
         self.distance = distance
         self.seed = seed
 
         if moontype == 'standard':
             N1, N2 = 1000, 1000
             data = utils.loadmat('pointclouds/two_moons')
-            Xin = data['features'][:dim].T
+            coords = data['features'][:dim].T
 
         elif moontype == 'synthesized':
             N1 = N // 2
             N2 = N - N1
 
-            coords1 = self._create_arc_moon(N1, sigmad, distance, 1, seed)
-            coords2 = self._create_arc_moon(N2, sigmad, distance, 2, seed)
+            coords1 = self._create_arc_moon(N1, variance, distance, 1, seed)
+            coords2 = self._create_arc_moon(N2, variance, distance, 2, seed)
 
-            Xin = np.concatenate((coords1, coords2))
+            coords = np.concatenate((coords1, coords2))
 
         else:
             raise ValueError('Unknown moontype {}'.format(moontype))
@@ -96,15 +96,14 @@ def __init__(self, moontype='standard', dim=2, sigmag=0.05,
             'vertex_size': 30,
         }
 
-        super(TwoMoons, self).__init__(Xin=Xin, sigma=sigmag, k=5,
-                                       center=False, rescale=False,
+        super(TwoMoons, self).__init__(coords, kernel_width=kernel_width, k=5,
                                        plotting=plotting, **kwargs)
 
     def _get_extra_repr(self):
         attrs = {'moontype': self.moontype,
                  'dim': self.dim,
-                 'sigmag': '{:.2f}'.format(self.sigmag),
-                 'sigmad': '{:.2f}'.format(self.sigmad),
+                 'kernel_width': '{:.2f}'.format(self.kernel_width),
+                 'variance': '{:.2f}'.format(self.variance),
                  'distance': '{:.2f}'.format(self.distance),
                  'seed': self.seed}
         attrs.update(super(TwoMoons, self)._get_extra_repr())
diff --git a/pygsp/reduction.py b/pygsp/reduction.py
index 50a8ae67..184d4b35 100644
--- a/pygsp/reduction.py
+++ b/pygsp/reduction.py
@@ -56,7 +56,7 @@ def graph_sparsify(M, epsilon, maxiter=10):
     Examples
     --------
     >>> from pygsp import reduction
-    >>> G = graphs.Sensor(256, k=20, distributed=True)
+    >>> G = graphs.Sensor(256, k=20, distributed=True, seed=0)
     >>> epsilon = 0.4
     >>> G2 = reduction.graph_sparsify(G, epsilon)
 
diff --git a/pygsp/tests/test_all.py b/pygsp/tests/test_all.py
index 2d52a185..1fc9a55d 100755
--- a/pygsp/tests/test_all.py
+++ b/pygsp/tests/test_all.py
@@ -14,9 +14,11 @@
 from pygsp.tests import test_docstrings
 from pygsp.tests import test_plotting
 from pygsp.tests import test_learning
+from pygsp.tests import test_nearest_neighbor
 
 
 suites = []
+suites.append(test_nearest_neighbor.suite)
 suites.append(test_graphs.suite)
 suites.append(test_filters.suite)
 suites.append(test_utils.suite)
diff --git a/pygsp/tests/test_filters.py b/pygsp/tests/test_filters.py
index 0bf066f5..cef1689d 100644
--- a/pygsp/tests/test_filters.py
+++ b/pygsp/tests/test_filters.py
@@ -149,7 +149,7 @@ def test_frame(self):
         def get_frame(freq_response):
             return self._G.U.dot(np.diag(freq_response).dot(self._G.U.T))
         gL = np.concatenate([get_frame(gl) for gl in g.evaluate(self._G.e)])
-        np.testing.assert_allclose(gL1, gL)
+        np.testing.assert_allclose(gL1, gL, atol=1e-10)
         np.testing.assert_allclose(gL2, gL, atol=1e-10)
 
     def test_complement(self, frame_bound=2.5):
diff --git a/pygsp/tests/test_graphs.py b/pygsp/tests/test_graphs.py
index d1e1c1af..63a9d7b9 100644
--- a/pygsp/tests/test_graphs.py
+++ b/pygsp/tests/test_graphs.py
@@ -9,6 +9,7 @@
 
 import os
 import random
+
 import sys
 import unittest
 
@@ -16,12 +17,15 @@
 import scipy.linalg
 from scipy import sparse
 import networkx as nx
-import graph_tool as gt
-import graph_tool.generation
 from skimage import data, img_as_float
 
 from pygsp import graphs
 
+import graph_tool as gt
+import graph_tool.generation
+
+
+
 
 class TestCase(unittest.TestCase):
 
@@ -448,6 +452,96 @@ def test_set_coordinates(self):
         G.set_coordinates('community2D')
         self.assertRaises(ValueError, G.set_coordinates, 'invalid')
 
+
+    def test_nngraph(self, n_vertices=24):
+        """Test all the combinations of metric, kind, backend."""
+        Graph = graphs.NNGraph
+        data = np.random.RandomState(42).uniform(size=(n_vertices, 3))
+        metrics = ['euclidean', 'manhattan', 'max_dist', 'minkowski']
+        # Not testing , 'flann', 'nmslib' as they are tested in test_nearest_neighbor
+        backends = ['scipy-kdtree', 'scipy-ckdtree'] 
+
+        for metric in metrics:
+            for kind in ['knn', 'radius']:
+                params = dict(features=data, metric=metric, kind=kind, k=6)
+                ref = Graph(backend='scipy-pdist', **params)
+                for backend in backends:
+                    graph = Graph(**params)
+                    np.testing.assert_allclose(graph.W.toarray(),
+                                                   ref.W.toarray(), rtol=1e-5)
+
+        # Distances between points on a circle.
+        angles = [0, 2 * np.pi / n_vertices]
+        points = np.stack([np.cos(angles), np.sin(angles)], axis=1)
+        distance = np.linalg.norm(points[0] - points[1])
+        weight = np.exp(-np.log(2) * distance**2)
+        column = np.zeros(n_vertices)
+        column[1] = weight
+        column[-1] = weight
+        adjacency = scipy.linalg.circulant(column)
+        data = graphs.Ring(n_vertices).coords
+        for kind in ['knn', 'radius']:
+            for backend in backends + ['scipy-pdist']:
+                graph = Graph(data, kind=kind, k=2, radius=1.01*distance,
+                              kernel_width=1, backend=backend)
+                np.testing.assert_allclose(graph.W.toarray(), adjacency)
+
+        graph = Graph(data, kind='radius', radius='estimate')
+        np.testing.assert_allclose(graph.radius, np.sqrt(8 / n_vertices))
+        graph = Graph(data, kind='radius', k=2, radius='estimate-knn')
+        np.testing.assert_allclose(graph.radius, distance)
+
+        graph = Graph(data, standardize=True)
+        np.testing.assert_allclose(np.mean(graph.coords, axis=0), 0, atol=1e-7)
+        np.testing.assert_allclose(np.std(graph.coords, axis=0), 1)
+
+        # Invalid parameters.
+        self.assertRaises(ValueError, Graph, np.ones(n_vertices))
+        self.assertRaises(ValueError, Graph, np.ones((n_vertices, 3, 4)))
+        self.assertRaises(ValueError, Graph, data, metric='invalid')
+        self.assertRaises(ValueError, Graph, data, kind='invalid')
+        self.assertRaises(ValueError, Graph, data, kernel='invalid')
+        self.assertRaises(ValueError, Graph, data, backend='invalid')
+        self.assertRaises(ValueError, Graph, data, kind='knn', k=0)
+        self.assertRaises(ValueError, Graph, data, kind='knn', k=n_vertices)
+        self.assertRaises(ValueError, Graph, data, kind='radius', radius=0)
+
+        # Empty graph.
+        if sys.version_info > (3, 4):  # no assertLogs in python 2.7
+            for backend in backends + ['scipy-pdist']:
+                if backend == 'nmslib':
+                    continue  # nmslib doesn't support radius
+                with self.assertLogs(level='WARNING'):
+                    graph = Graph(data, kind='radius', radius=1e-9,
+                                  backend=backend)
+                self.assertEqual(graph.n_edges, 0)
+
+        # Backend parameters.
+        Graph(data, lap_type='normalized')
+        Graph(data, plotting=dict(vertex_size=10))
+        Graph(data, backend='flann', algorithm='kmeans')
+        Graph(data, backend='nmslib', method='vptree')
+        Graph(data, backend='nmslib', index=dict(post=2))
+        Graph(data, backend='nmslib', query=dict(efSearch=100))
+        for backend in ['scipy-kdtree', 'scipy-ckdtree']:
+            Graph(data, backend=backend, eps=1e-2)
+            Graph(data, backend=backend, leafsize=9)
+        self.assertRaises(ValueError, Graph, data, backend='scipy-pdist', a=0)
+
+        # Kernels.
+        for name, kernel in graphs.NNGraph._kernels.items():
+            similarity = 0 if name == 'rectangular' else 0.5
+            np.testing.assert_allclose(kernel(np.ones(10)), similarity)
+            np.testing.assert_allclose(kernel(np.zeros(10)), 1)
+        Graph(data, kernel=lambda d: d.min()/d)
+        if sys.version_info > (3, 4):  # no assertLogs in python 2.7
+            with self.assertLogs(level='WARNING'):
+                Graph(data, kernel=lambda d: 1/d)
+
+        # Attributes.
+        self.assertEqual(Graph(data, kind='knn').radius, None)
+        self.assertEqual(Graph(data, kind='radius').k, None)
+
     def test_line_graph(self):
         adjacency = [
             [0, 1, 1, 3],
@@ -492,24 +586,6 @@ def test_subgraph(self, n_vertices=100):
         self.assertIs(graph.lap_type, self._G.lap_type)
         self.assertEqual(graph.plotting, self._G.plotting)
 
-    def test_nngraph(self, n_vertices=30):
-        rs = np.random.RandomState(42)
-        Xin = rs.normal(size=(n_vertices, 3))
-        dist_types = ['euclidean', 'manhattan', 'max_dist', 'minkowski']
-
-        for dist_type in dist_types:
-
-            # Only p-norms with 1<=p<=infinity permitted.
-            if dist_type != 'minkowski':
-                graphs.NNGraph(Xin, NNtype='radius', dist_type=dist_type)
-                graphs.NNGraph(Xin, NNtype='knn', dist_type=dist_type)
-
-            # Distance type unsupported in the C bindings,
-            # use the C++ bindings instead.
-            if dist_type != 'max_dist':
-                graphs.NNGraph(Xin, use_flann=True, NNtype='knn',
-                               dist_type=dist_type)
-
     def test_bunny(self):
         graphs.Bunny()
 
diff --git a/pygsp/tests/test_nearest_neighbor.py b/pygsp/tests/test_nearest_neighbor.py
new file mode 100644
index 00000000..c11d38b7
--- /dev/null
+++ b/pygsp/tests/test_nearest_neighbor.py
@@ -0,0 +1,79 @@
+import unittest
+import numpy as np
+from pygsp._nearest_neighbor import nearest_neighbor, sparse_distance_matrix
+
+class TestCase(unittest.TestCase):
+    def test_nngraph(self, n_vertices=100):
+        data = np.random.RandomState(42).uniform(size=(n_vertices, 3))
+        metrics = ['euclidean', 'manhattan', 'max_dist', 'minkowski']
+        backends = ['scipy-kdtree', 'scipy-ckdtree', 'flann', 'nmslib']
+
+        for metric in metrics:
+            for kind in ['knn', 'radius']:
+                for backend in backends:
+                    params = dict(features=data, metric=metric, kind=kind, radius=0.25, k=8)
+                    
+                    ref_nn, ref_d = nearest_neighbor(backend='scipy-pdist', **params)
+                    # Unsupported combinations.
+                    if backend == 'flann' and metric == 'max_dist':
+                        self.assertRaises(ValueError, nearest_neighbor, data,
+                                          metric=metric, backend=backend)
+                    elif backend == 'nmslib' and metric == 'minkowski':
+                        self.assertRaises(ValueError, nearest_neighbor, data,
+                                          metric=metric, backend=backend)
+                    elif backend == 'nmslib' and kind == 'radius':
+                        self.assertRaises(ValueError, nearest_neighbor, data,
+                                          kind=kind, backend=backend)
+                    else:
+                        params['backend'] = backend
+                        if backend == 'flann':
+                            other_nn, other_d = nearest_neighbor(random_seed=44, target_precision=1, **params)
+                        else:
+                            other_nn, other_d = nearest_neighbor(**params)
+                        print(kind, backend)
+                        if backend == 'flann':
+                            for a,b in zip(ref_nn, other_nn):
+                                assert(len(set(a)-set(b))<=2)
+
+                            for a,b in zip(ref_d, other_d):
+                                np.testing.assert_allclose(np.mean(a).astype(np.float32),np.mean(b), atol=2e-2)                            
+                        else:
+                            for a,b,c,d in zip(ref_nn, other_nn, ref_d, other_d):
+                                e = set(a)-set(b)
+                                assert(len(e)<=1)
+                                if len(e)==0:
+                                    np.testing.assert_allclose(np.sort(c),np.sort(d), rtol=1e-5)
+
+    def test_sparse_distance_matrix(self):
+        data = np.random.RandomState(42).uniform(size=(200, 3))
+        neighbors, distances = nearest_neighbor(data)
+        W = sparse_distance_matrix(neighbors, distances, symmetrize=True)
+        # Assert symetry
+        np.testing.assert_allclose(W.todense(), W.T.todense())
+        # positivity
+        np.testing.assert_array_equal(W.todense()>=0, True)
+        # 0 diag
+        np.testing.assert_array_equal(np.diag(W.todense())==0, True)
+
+        # Assert that it is not symmetric anymore
+        W = sparse_distance_matrix(neighbors, distances, symmetrize=False)
+        assert(np.sum(np.abs(W.todense()-W.T.todense()))>0.1)
+        # positivity
+        np.testing.assert_array_equal(W.todense()>=0, True)
+        # 0 diag
+        np.testing.assert_array_equal(np.diag(W.todense())==0, True)
+        # everything is used once
+        np.testing.assert_allclose(np.sum(W.todense()), np.sum(distances))
+
+        # simple test with a kernel
+        W = sparse_distance_matrix(neighbors, 1/(1+distances), symmetrize=True)
+        # Assert symetry
+        np.testing.assert_allclose(W.todense(), W.T.todense())
+        # positivity
+        np.testing.assert_array_equal(W.todense()>=0, True)
+        # 0 diag
+        np.testing.assert_array_equal(np.diag(W.todense())==0, True)
+
+
+suite = unittest.TestLoader().loadTestsFromTestCase(TestCase)
+
diff --git a/pygsp/tests/test_plotting.py b/pygsp/tests/test_plotting.py
index 822d3127..3694b4f7 100644
--- a/pygsp/tests/test_plotting.py
+++ b/pygsp/tests/test_plotting.py
@@ -50,8 +50,8 @@ def test_plot_graphs(self):
 
             # Classes who require parameters.
             if classname == 'NNGraph':
-                Xin = np.arange(90).reshape(30, 3)
-                Gs.append(Graph(Xin))
+                features = np.random.RandomState(42).normal(size=(30, 3))
+                Gs.append(Graph(features))
             elif classname in ['ImgPatches', 'Grid2dImgPatches']:
                 Gs.append(Graph(img=self._img, patch_shape=(3, 3)))
             elif classname == 'LineGraph':
diff --git a/setup.py b/setup.py
index 2ed970f9..37849940 100644
--- a/setup.py
+++ b/setup.py
@@ -40,8 +40,8 @@
             # Construct patch graphs from images.
             'scikit-image',
             # Approximate nearest neighbors for kNN graphs.
-            'pyflann; python_version == "2.*"',
-            'pyflann3; python_version == "3.*"',
+            'cyflann',
+            'nmslib',
             # Convex optimization on graph.
             'pyunlocbox',
             # Plot graphs, signals, and filters.