.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_examples/plot_03_spectrum.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. .. rst-class:: sphx-glr-example-title .. _sphx_glr_auto_examples_plot_03_spectrum.py: The Laplacian and its spectrum ============================== The Laplace-Beltrami operator is the single object the rest of pyFM is built on. Its eigenfunctions form an orthonormal basis of the functions living on a surface, ordered from the slowest to the fastest varying. Truncating that basis is what turns a correspondence between two shapes into a small matrix called a functional map. This page builds the operator, looks at its spectrum, differentiates and integrates with it, compresses functions in its basis, and ends on the descriptors that the functional-map examples use as input. .. GENERATED FROM PYTHON SOURCE LINES 16-20 Locating data ------------- The example meshes are stored in ``examples/data`` and are not part of the package. .. GENERATED FROM PYTHON SOURCE LINES 20-48 .. code-block:: Python from pathlib import Path import matplotlib.pyplot as plt import numpy as np import pyvista as pv from scipy import sparse import pyFM.signatures as sg from pyFM.mesh import TriMesh from pyFM.viz import plot_mesh def example_data(name): """Return the path to a mesh in the ``examples/data`` folder.""" for folder in [Path.cwd(), *Path.cwd().parents]: candidate = folder / "examples" / "data" / name if candidate.exists(): return candidate raise FileNotFoundError( f"{name!r} not found. These examples read data from the pyFM repository; " "clone it and run them there." ) mesh = TriMesh.load(example_data("cat.obj"), area_normalize=True, center=True) .. GENERATED FROM PYTHON SOURCE LINES 50-69 Building the operator --------------------- :meth:`~pyFM.mesh.trimesh.TriMesh.process` discretizes the Laplacian and computes the first ``k`` eigenpairs, storing everything on the mesh object. It returns the mesh itself, so it chains, and it never recomputes a spectrum it already has (asking for fewer eigenvectors than are stored simply truncates it). The discrete Laplacian is usually defined as :math:`\Delta = A^{-1} W`, where :math:`W` is the stiffness matrix and :math:`A` is the mass matrix. The eigenproblem is therefore :math:`W \phi_k = \lambda_k A \phi_k`, which is a generalized eigenproblem. The eigenvectors are orthonormal with respect to the mass matrix, :math:`\phi^T A \phi_j = \delta_{ij}`. The default discretization is the standard cotangent scheme. ``intrinsic=True`` switches to an intrinsic triangulation, which ensures positive cotangent weights. ``robust=True`` switches to the tufted Laplacian, which tolerates non-manifold input and works for point clouds. Using ``k=0``, or calling :meth:`~pyFM.mesh.trimesh.TriMesh.compute_operators`, builds the matrices without solving the eigenproblem. .. GENERATED FROM PYTHON SOURCE LINES 69-77 .. code-block:: Python mesh.process(k=100) print(f"stiffness: {mesh.stiffness.shape}, {mesh.stiffness.nnz} nonzeros") print(f"mass : {mesh.mass.shape}, {mesh.mass.nnz} nonzeros") print(f"eigenvalues : {mesh.eigenvalues.shape}") print(f"eigenvectors: {mesh.eigenvectors.shape}") .. rst-class:: sphx-glr-script-out .. code-block:: none stiffness: (7207, 7207), 50437 nonzeros mass : (7207, 7207), 7207 nonzeros eigenvalues : (100,) eigenvectors: (7207, 100) .. GENERATED FROM PYTHON SOURCE LINES 78-81 The two matrices are available under their full names, ``stiffness`` and ``mass``, and under the short aliases ``W`` and ``A`` that keep formulas simple. The mass matrix is the discrete integration weight, so summing it against a constant function recovers the area of the shape. .. GENERATED FROM PYTHON SOURCE LINES 81-88 .. code-block:: Python ones = np.ones(mesh.n_vertices) print(f"mesh.W is mesh.stiffness: {mesh.W is mesh.stiffness}") print(f"integral of 1 : {mesh.integrate(ones):.6f}") print(f"mesh area : {mesh.area:.6f}") .. rst-class:: sphx-glr-script-out .. code-block:: none mesh.W is mesh.stiffness: True integral of 1 : 1.000000 mesh area : 1.000000 .. GENERATED FROM PYTHON SOURCE LINES 89-92 The eigenvectors are orthonormal *with respect to the mass matrix*, not the plain dot product. That will make :meth:`~pyFM.mesh.trimesh.TriMesh.project` below a simple matrix product and not a linear solve. .. GENERATED FROM PYTHON SOURCE LINES 92-97 .. code-block:: Python gram = mesh.eigenvectors.T @ mesh.mass @ mesh.eigenvectors print(f"max |V^T A V - I| = {np.abs(gram - np.eye(gram.shape[0])).max():.2e}") .. rst-class:: sphx-glr-script-out .. code-block:: none max |V^T A V - I| = 4.44e-15 .. GENERATED FROM PYTHON SOURCE LINES 98-105 The spectrum ------------ The eigenvalues are non-negative and sorted. The first one is zero, with the constant function as eigenvector. Then an eigenvalue is a frequency: the larger it is, the faster its eigenfunction oscillates. This is a generalization of the Fourier basis to curved surfaces, and it is the reason why the span of the first eigenfunctions is used to compress signals. .. GENERATED FROM PYTHON SOURCE LINES 105-117 .. code-block:: Python print(f"lambda_0 = {mesh.eigenvalues[0]:.2e} (zero, up to numerical error)") print(f"next eigenvalues: {np.round(mesh.eigenvalues[1:5], 4)}") fig, ax = plt.subplots(figsize=(6, 3.5), constrained_layout=True) ax.plot(mesh.eigenvalues, ".-", markersize=4) ax.set_xlabel("index $k$") ax.set_ylabel(r"eigenvalue $\lambda_k$") ax.set_title("Laplace-Beltrami spectrum") ax.grid(alpha=0.3) plt.show() .. image-sg:: /auto_examples/images/sphx_glr_plot_03_spectrum_001.png :alt: Laplace-Beltrami spectrum :srcset: /auto_examples/images/sphx_glr_plot_03_spectrum_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none lambda_0 = 6.30e-14 (zero, up to numerical error) next eigenvalues: [ 6.2723 11.9226 18.58 23.4044] .. GENERATED FROM PYTHON SOURCE LINES 118-123 Eigenfunctions -------------- Each eigenvector is a function on the vertices, plotted here with a diverging colormap. Note that higher frequency eigenvectors split the mesh into smaller regions. .. GENERATED FROM PYTHON SOURCE LINES 123-136 .. code-block:: Python ev_inds = [1, 5, 20] pl = pv.Plotter(shape=(1, len(ev_inds))) for i, ev_ind in enumerate(ev_inds): pl.subplot(0, i) evec = mesh.eigenvectors[:, ev_ind] bound = np.abs(evec).max() plot_mesh(mesh, scalars=evec, cmap="coolwarm", clim=(-bound, bound), pl=pl) pl.add_title(f"mode {ev_ind}, lambda = {mesh.eigenvalues[ev_ind]:.1f}", font_size=9) pl.link_views() pl.show() .. tab-set:: .. tab-item:: Static Scene .. image-sg:: /auto_examples/images/sphx_glr_plot_03_spectrum_002.png :alt: plot 03 spectrum :srcset: /auto_examples/images/sphx_glr_plot_03_spectrum_002.png :class: sphx-glr-single-img .. tab-item:: Interactive Scene .. offlineviewer:: /home/runner/work/pyFM/pyFM/docs/source/auto_examples/images/sphx_glr_plot_03_spectrum_002.vtksz .. GENERATED FROM PYTHON SOURCE LINES 137-144 Gradient and divergence ----------------------- :meth:`~pyFM.mesh.trimesh.TriMesh.gradient` takes a function on the ``n`` vertices and returns a vector field on the ``m`` faces: on each triangle the function is linear, so its gradient is a single vector tangent to that triangle. A constant function therefore has a strictly zero gradient. .. GENERATED FROM PYTHON SOURCE LINES 144-164 .. code-block:: Python f = mesh.eigenvectors[:, 1] grad_f = mesh.gradient(f) print(f"function: {f.shape} (per vertex) -> gradient: {grad_f.shape} (per face)") print(f"max |gradient of a constant|: {np.abs(mesh.gradient(ones)).max():.2e}") bound = np.abs(f).max() plot_mesh( mesh, scalars=f, cmap="coolwarm", clim=(-bound, bound), vfield=grad_f, vfield_rescale=0.01, vfield_tolerance=0.01, vfield_color="black", ) .. tab-set:: .. tab-item:: Static Scene .. image-sg:: /auto_examples/images/sphx_glr_plot_03_spectrum_003.png :alt: plot 03 spectrum :srcset: /auto_examples/images/sphx_glr_plot_03_spectrum_003.png :class: sphx-glr-single-img .. tab-item:: Interactive Scene .. offlineviewer:: /home/runner/work/pyFM/pyFM/docs/source/auto_examples/images/sphx_glr_plot_03_spectrum_003.vtksz .. rst-class:: sphx-glr-script-out .. code-block:: none function: (7207,) (per vertex) -> gradient: (14410, 3) (per face) max |gradient of a constant|: 0.00e+00 .. GENERATED FROM PYTHON SOURCE LINES 165-175 :meth:`~pyFM.mesh.trimesh.TriMesh.divergence` goes from a per-face vector field back to a function on the vertices. Composing the two gives the Laplacian itself, :math:`\Delta = \text{div} \circ \text{grad}`. We can check this numerically, by verifying :math:`\text{div}(\text{grad} \, \phi_k) = \lambda_k \phi_k` for an eigenfunction :math:`\phi_k`. .. warning:: pyFM uses the *geometry processing* sign convention, where the Laplacian :math:`\Delta = A^{-1} W` has a non-negative spectrum. ``divergence(gradient(f))`` therefore returns :math:`+\Delta f`, the opposite sign to the usual math convention. .. GENERATED FROM PYTHON SOURCE LINES 175-185 .. code-block:: Python div_grad_f = mesh.divergence(grad_f) lambda_1 = mesh.eigenvalues[1] rel_error = np.linalg.norm(div_grad_f - lambda_1 * f) / np.linalg.norm(lambda_1 * f) print(f"eigenvalue lambda_1 : {lambda_1:.10f}") print(f"median of div(grad f) / f : {np.median(div_grad_f / f):.10f}") print(f"relative error on Delta f - lambda f: {rel_error:.2e}") .. rst-class:: sphx-glr-script-out .. code-block:: none eigenvalue lambda_1 : 6.2723447940 median of div(grad f) / f : 6.2723447940 relative error on Delta f - lambda f: 2.48e-10 .. GENERATED FROM PYTHON SOURCE LINES 186-197 Compressing functions in the spectral basis ------------------------------------------- This is the idea the whole functional map framework rests on. Because the eigenfunctions form an orthonormal basis, :meth:`~pyFM.mesh.trimesh.TriMesh.project` turns a function over the ``n`` vertices into a few coefficients. :meth:`~pyFM.mesh.trimesh.TriMesh.unproject` turns them back into values on the surface, and :meth:`~pyFM.mesh.trimesh.TriMesh.reconstruct` does both at once. We generate a function on the surface by transforming a geodesic distance field. .. GENERATED FROM PYTHON SOURCE LINES 197-204 .. code-block:: Python target = np.cos(2 * np.pi * mesh.geodesic_from(200) / 0.5) coefficients = mesh.project(target) print(f"function: {target.shape} -> coefficients: {coefficients.shape}") print(f"first coefficients: {np.round(coefficients[:5], 4)}") .. rst-class:: sphx-glr-script-out .. code-block:: none function: (7207,) -> coefficients: (100,) first coefficients: [-0.1536 -0.0226 -0.0622 0.0329 -0.0046] .. GENERATED FROM PYTHON SOURCE LINES 205-208 Projection turns a function into K coefficients. Higher values of K allow for a more accurate reconstruction, but at the cost of storing more data. We see that this field is hard to reproduce when using too few eigenvectors. However, with 100 coefficients, we are able to reconstruct the function with good accuracy, independently from the number of vertices. .. GENERATED FROM PYTHON SOURCE LINES 208-223 .. code-block:: Python k_values = [5, 20, 100] clim = (target.min(), target.max()) pl = pv.Plotter(shape=(1, 1 + len(k_values))) pl.subplot(0, 0) plot_mesh(mesh, scalars=target, cmap="plasma", clim=clim, pl=pl) pl.add_title("original", font_size=9) for i, k in enumerate(k_values): pl.subplot(0, i + 1) plot_mesh(mesh, scalars=mesh.reconstruct(target, k=k), cmap="plasma", clim=clim, pl=pl) pl.add_title(f"k = {k}", font_size=9) pl.link_views() pl.show() .. tab-set:: .. tab-item:: Static Scene .. image-sg:: /auto_examples/images/sphx_glr_plot_03_spectrum_004.png :alt: plot 03 spectrum :srcset: /auto_examples/images/sphx_glr_plot_03_spectrum_004.png :class: sphx-glr-single-img .. tab-item:: Interactive Scene .. offlineviewer:: /home/runner/work/pyFM/pyFM/docs/source/auto_examples/images/sphx_glr_plot_03_spectrum_004.vtksz .. GENERATED FROM PYTHON SOURCE LINES 224-225 :meth:`~pyFM.mesh.trimesh.TriMesh.l2_sqnorm` measures the error by integrating on the surface (weighting with vertex area) .. GENERATED FROM PYTHON SOURCE LINES 225-242 .. code-block:: Python k_range = np.arange(1, mesh.eigenvalues.size + 1) errors = [ np.sqrt(mesh.l2_sqnorm(target - mesh.reconstruct(target, k=k)) / mesh.l2_sqnorm(target)) for k in k_range ] fig, ax = plt.subplots(figsize=(6, 3.5), constrained_layout=True) ax.semilogy(k_range, errors, ".-", markersize=4) ax.set_xlabel("number of eigenfunctions $k$") ax.set_ylabel("relative $L^2$ error") ax.set_title("Reconstruction of a geodesic distance field") ax.grid(alpha=0.3, which="both") plt.show() print(f"relative L2 error with k = 100: {errors[-1]:.3f}") .. image-sg:: /auto_examples/images/sphx_glr_plot_03_spectrum_005.png :alt: Reconstruction of a geodesic distance field :srcset: /auto_examples/images/sphx_glr_plot_03_spectrum_005.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none relative L2 error with k = 100: 0.083 .. GENERATED FROM PYTHON SOURCE LINES 243-264 Diffusing heat -------------- Heat spreading over the surface obeys .. math:: \frac{\partial u}{\partial t} = - \Delta u The minus sign goes with the convention that: :math:`\Delta = A^{-1} W` is positive semi-definite. In the spectral basis the equation decouples, and each coefficient simply decays as :math:`e^{-\lambda_k t}`. Discretizing the time derivative explicitly (forward Euler) gives .. math:: u_{n+1} = u_n - \delta t \, A^{-1} W u_n One sparse product per step, but stable only while :math:`\delta t < 2 / \lambda_{\max}`. Since the *largest* eigenvalue scales like the inverse squared edge length, this is an issue for badly triangulated shapes, which is why implicit methods are often preferred. .. GENERATED FROM PYTHON SOURCE LINES 264-272 .. code-block:: Python lambda_max = sparse.linalg.eigsh( mesh.W.tocsc(), k=1, M=mesh.A.tocsc(), which="LM", return_eigenvectors=False )[0] print(f"largest eigenvalue: {lambda_max:.2e}") print(f"explicit steps must satisfy dt < {2 / lambda_max:.2e}") .. rst-class:: sphx-glr-script-out .. code-block:: none largest eigenvalue: 3.99e+07 explicit steps must satisfy dt < 5.01e-08 .. GENERATED FROM PYTHON SOURCE LINES 273-279 The implicit (backward Euler) step instead solves .. math:: (A + \delta t \, W) \, u_{n+1} = A \, u_n which is unconditionally stable. .. GENERATED FROM PYTHON SOURCE LINES 279-290 .. code-block:: Python source = 200 delta = np.zeros(mesh.n_vertices) delta[source] = 1.0 # this is A @ dirac; the area weight cancels out diffusion_time = 1e0 heat = sparse.linalg.spsolve(mesh.A + diffusion_time * mesh.W, delta) print(f"diffusion time: {diffusion_time:.0e}, one linear solve") print(f"integral of the solution: {mesh.integrate(heat):.4f} (heat is conserved)") .. rst-class:: sphx-glr-script-out .. code-block:: none diffusion time: 1e+00, one linear solve integral of the solution: 1.0000 (heat is conserved) .. GENERATED FROM PYTHON SOURCE LINES 291-292 We plot the solution .. GENERATED FROM PYTHON SOURCE LINES 292-301 .. code-block:: Python plot_mesh( mesh, scalars=heat, cmap="Reds", points=source, points_color="black", ) .. tab-set:: .. tab-item:: Static Scene .. image-sg:: /auto_examples/images/sphx_glr_plot_03_spectrum_006.png :alt: plot 03 spectrum :srcset: /auto_examples/images/sphx_glr_plot_03_spectrum_006.png :class: sphx-glr-single-img .. tab-item:: Interactive Scene .. offlineviewer:: /home/runner/work/pyFM/pyFM/docs/source/auto_examples/images/sphx_glr_plot_03_spectrum_006.vtksz .. GENERATED FROM PYTHON SOURCE LINES 302-313 Spectral descriptors -------------------- The Heat and Wave Kernel Signatures are built directly from the spectrum. Both give every vertex a vector of values that describes the geometry around it, and both are *isometry invariant*: they depend on the surface, not on how it is posed in space. That is what makes them usable to put two different shapes in correspondence. :func:`~pyFM.signatures.HKS_functions.mesh_HKS` and :func:`~pyFM.signatures.WKS_functions.mesh_WKS` read the spectrum off an already processed mesh and pick their own time and energy scales from it. .. GENERATED FROM PYTHON SOURCE LINES 313-319 .. code-block:: Python hks = sg.mesh_HKS(mesh, num_T=100) wks = sg.mesh_WKS(mesh, num_E=100) print(f"HKS: {hks.shape}, WKS: {wks.shape}") .. rst-class:: sphx-glr-script-out .. code-block:: none HKS: (7207, 100), WKS: (7207, 100) .. GENERATED FROM PYTHON SOURCE LINES 320-322 The HKS describes how much heat remains at a vertex after diffusing for a time ``t``. Short times see only the local curvature, while longer times see larger scales. .. GENERATED FROM PYTHON SOURCE LINES 322-333 .. code-block:: Python time_inds = [0, 30, 80] pl = pv.Plotter(shape=(1, len(time_inds))) for i, t_ind in enumerate(time_inds): pl.subplot(0, i) plot_mesh(mesh, scalars=hks[:, t_ind], cmap="coolwarm", pl=pl) pl.add_title(f"HKS, time index {t_ind}", font_size=9) pl.link_views() pl.show() .. tab-set:: .. tab-item:: Static Scene .. image-sg:: /auto_examples/images/sphx_glr_plot_03_spectrum_007.png :alt: plot 03 spectrum :srcset: /auto_examples/images/sphx_glr_plot_03_spectrum_007.png :class: sphx-glr-single-img .. tab-item:: Interactive Scene .. offlineviewer:: /home/runner/work/pyFM/pyFM/docs/source/auto_examples/images/sphx_glr_plot_03_spectrum_007.vtksz .. GENERATED FROM PYTHON SOURCE LINES 334-336 The WKS replaces the diffusion time by an energy band, which makes it more selective: each column responds to one range of frequencies. .. GENERATED FROM PYTHON SOURCE LINES 336-347 .. code-block:: Python energy_inds = [0, 30, 80] pl = pv.Plotter(shape=(1, len(energy_inds))) for i, e_ind in enumerate(energy_inds): pl.subplot(0, i) plot_mesh(mesh, scalars=wks[:, e_ind], cmap="coolwarm", pl=pl) pl.add_title(f"WKS, energy index {e_ind}", font_size=9) pl.link_views() pl.show() .. tab-set:: .. tab-item:: Static Scene .. image-sg:: /auto_examples/images/sphx_glr_plot_03_spectrum_008.png :alt: plot 03 spectrum :srcset: /auto_examples/images/sphx_glr_plot_03_spectrum_008.png :class: sphx-glr-single-img .. tab-item:: Interactive Scene .. offlineviewer:: /home/runner/work/pyFM/pyFM/docs/source/auto_examples/images/sphx_glr_plot_03_spectrum_008.vtksz .. GENERATED FROM PYTHON SOURCE LINES 348-352 These are exactly the arrays :class:`~pyFM.functional.FunctionalMapping` computes when given ``descr_type="HKS"`` or ``descr_type="WKS"``: a functional map is optimized so that it maps the descriptors of one shape onto the descriptors of the other, expressed in the truncated spectral basis built above. .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 3.180 seconds) .. _sphx_glr_download_auto_examples_plot_03_spectrum.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot_03_spectrum.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_03_spectrum.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: plot_03_spectrum.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_