diff --git a/geomfum/descriptor/spatial.py b/geomfum/descriptor/spatial.py
new file mode 100644
index 00000000..8a7aab6b
--- /dev/null
+++ b/geomfum/descriptor/spatial.py
@@ -0,0 +1,73 @@
+import numpy as np
+import scipy
+from sklearn.preprocessing import normalize
+
+from geomfum.sample import NearestNeighborsIndexSampler
+
+
+def linear_compact(x):
+    """Linearly decreasing function between 0 and 1."""
+    return 1 - x
+
+
+def poly_compact(x):
+    """Polynomial decreasing function between 0 and 1."""
+    return 1 - 3 * x**2 + 2 * x**3
+
+
+def exp_compact(x):
+    """Exponential decreasing function between 0 and 1."""
+    return np.exp(1 - 1 / (1 - x**2))
+
+
+class NasikunLocalFunctionsConstructor:
+    """
+
+    https://onlinelibrary.wiley.com/doi/full/10.1111/cgf.13496
+
+    Parameters
+    ----------
+    min_n_samples : int
+        Minimum number of samples to target.
+        Ignored if ``sampler is not None``.
+    """
+
+    # TODO: think about inheritance
+    # TODO: check if this is the original source
+    # TODO: still a naive implementation
+
+    def __init__(self, min_n_samples=100, sampler=None, local_transform=None):
+        if sampler is None:
+            sampler = NearestNeighborsIndexSampler(min_n_samples=min_n_samples)
+
+        if local_transform is None:
+            local_transform = poly_compact
+
+        self.sampler = sampler
+        self.local_transform = local_transform
+
+    def __call__(self, shape):
+        # assumes equipped shape
+
+        # TODO: think about how to pass metric
+        # TODO: need to implement clone for metric, in case of adaptation
+        subsample = self.sampler.sample(shape)
+
+        dists, targets = shape.metric.dist(subsample)
+
+        data = np.hstack(dists)
+        row_indices = []
+        column_indices = []
+        for index, targets_ in enumerate(targets):
+            row_indices.extend(targets_.tolist())
+            column_indices.extend([index] * len(targets_))
+
+        matrix = scipy.sparse.csr_matrix(
+            (data, (row_indices, column_indices)),
+            shape=(shape.n_vertices, subsample.size),
+        )
+        matrix.data = self.local_transform(matrix.data)
+
+        matrix = normalize(matrix, axis=1, norm="l1")
+
+        return matrix
diff --git a/geomfum/laplacian.py b/geomfum/laplacian.py
index 6a05bf8a..fee82b8b 100644
--- a/geomfum/laplacian.py
+++ b/geomfum/laplacian.py
@@ -5,6 +5,7 @@
 import geomfum.wrap as _wrap  # noqa (for register)
 from geomfum._registry import LaplacianFinderRegistry, MeshWhichRegistryMixins
 from geomfum.basis import LaplaceEigenBasis
+from geomfum.descriptor.spatial import NasikunLocalFunctionsConstructor
 from geomfum.numerics.eig import ScipyEigsh
 
 
@@ -131,3 +132,104 @@ def __call__(self, shape, as_basis=True, recompute=False):
             return LaplaceEigenBasis(shape, eigenvals, eigenvecs)
 
         return eigenvals, eigenvecs
+
+
+class NasikunLaplacianSpectrumFinder:
+    """Algorithm to find Laplacian spectrum.
+
+    https://onlinelibrary.wiley.com/doi/full/10.1111/cgf.13496
+
+    Parameters
+    ----------
+    spectrum_size : int
+        Spectrum size. Ignored if ``eig_solver`` is not None.
+    nonzero : bool
+        Remove zero zero eigenvalue.
+    fix_sign : bool
+        Wheather to have all the first components with positive sign.
+    laplacian_finder : BaseLaplacianFinder
+        Algorithm to find the Laplacian. Ignored if Laplace and mass matrices
+        were already computed.
+    eig_solver : EigSolver
+        Eigen solver.
+    """
+
+    # TODO: can be used under the hierarchical case
+    # TODO: update inheritance structure
+
+    def __init__(
+        self,
+        spectrum_size=100,
+        nonzero=False,
+        fix_sign=False,
+        laplacian_finder=None,
+        eig_solver=None,
+        local_func_constr=None,
+        min_n_samples=150,
+    ):
+        # min_n_samples ignored if local_func_constr is not None
+        if eig_solver is None:
+            eig_solver = ScipyEigsh(spectrum_size=spectrum_size, sigma=-0.01)
+
+        if local_func_constr is None:
+            local_func_constr = NasikunLocalFunctionsConstructor(
+                min_n_samples=min_n_samples
+            )
+
+        self.nonzero = nonzero
+        self.fix_sign = fix_sign
+        self.laplacian_finder = laplacian_finder
+        self.eig_solver = eig_solver
+        self.local_func_constr = local_func_constr
+
+    def __call__(self, shape, as_basis=True, recompute=False):
+        """Apply algorithm.
+
+        Parameters
+        ----------
+        shape : Shape
+            Shape.
+        as_basis : bool
+            Whether return basis or eigenvals/vecs.
+        recompute : bool
+            Whether to recompute Laplacian if information is cached.
+
+        Returns
+        -------
+        eigenvals : array-like, shape=[spectrum_size]
+            Eigenvalues. (If ``basis is False``.)
+        eigenvecs : array-like, shape=[n_vertices, spectrum_size]
+            Eigenvectors. (If ``basis is False``.)
+        basis : LaplaceEigenBasis
+            A basis. (If ``basis is True``.)
+        """
+        stiffness_matrix, mass_matrix = shape.laplacian.find(
+            self.laplacian_finder, recompute=recompute
+        )
+
+        local_func_mat = self.local_func_constr(shape)
+        restricted_mass_matrix = (
+            local_func_mat.T @ shape.laplacian.mass_matrix @ local_func_mat
+        )
+        restricted_stiffness_matrix = (
+            local_func_mat.T @ shape.laplacian.stiffness_matrix @ local_func_mat
+        )
+
+        eigenvals, restricted_eigenvecs = self.eig_solver(
+            restricted_stiffness_matrix, M=restricted_mass_matrix
+        )
+        eigenvecs = local_func_mat @ restricted_eigenvecs
+
+        if self.nonzero:
+            eigenvals = eigenvals[1:]
+            eigenvecs = eigenvecs[:, 1:]
+
+        if self.fix_sign:
+            indices = eigenvecs[0, :] < 0
+            eigenvals[indices] *= -1
+            eigenvecs[:, indices] *= -1
+
+        if as_basis:
+            return LaplaceEigenBasis(shape, eigenvals, eigenvecs)
+
+        return eigenvals, eigenvecs