Skip to content

ManifoldHarmonics

Jump to bottom
Bruno Levy edited this page Mar 13, 2022 · 4 revisions

Manifold Harmonics

Introduction

ManifoldHarmonics defines a Fourier-like basis on arbitrary meshes [1],[2]. Let us consider the classical Fourier basis. It decomposes a signal on a sum of sines and cosines. One can think about the sine and cosine function as the eigenfunction of a linear operator. Let us see what it means:

  • an operator is a "function of functions", that is, something that takes as an argument a function, and that returns another function. If an operator is denoted by L, its application to a function f is denoted by Lf (and it is another function). For instance, derivation is an operator;

  • an operator L is said to be linear if L(af+bg) = aLf + bLg for any couple of functions f and g, and any couple of scalars a and b. For instance, derivation is a linear operator;

  • given a linear operator L, a function f is said to be an eigenfunction of L if there exists a scalar a such that Lf = af, that is, applying L to f scales f by a factor a. The scalar a is called an eigenvalue of L. For instance, the function f defined by f(x) = exp(ax) is an eigenfunction of derivation, and the associated eigenvalue is a.

Let us now consider second order derivative and the function f defined by f(x) = cos(ax+b). The second order derivative of f is defined by f''(x) = -(a^2)cos(ax+b), hence f is an eigenfunction of the second order derivative associated with the eigenvalue -(a^2). This gives us a particular point of view of the sines and cosines functions used by the Fourier transform. They are the eigenfunctions of the second-order derivative. Moreover, it is easy to check that they form an orthogonal basis with respect to a certain "dot product" for functions (see [1],[2] for more details).

With this particular point of view, we can "port" this idea to more general settings, by "porting" the notion of second-order derivatives to more general settings:

  • in 2D, it is the Laplace operator, that is the sum of the second-order derivatives with respect to x and to y, that plays the role of the second-order derivative. The eigenfunctions of the Laplace operator define the Discrete Cosine Transform, widely used in image compression. They also correspond to stationnary waves on a vibrating plate, characterized by Ernst Chladni.

  • Another formula for the Laplace operator is given by Lf = div(grad(f)). This formula is interesting, because it involves the gradient and the divergence, that can be easily generalized to curved surfaces (Laplace-Beltrami operator). Doing so on the sphere, one can compute the eigenfunctions in closed form. This defines spherical harmonics.

  • On more general surfaces, one can compute the eigenfunctions of the Laplace-Beltrami operator, using the Finite Elements approach (Manifold Harmonics). This boils down to compute the eigenvectors of a matrix, which can be implemented using numerical methods (such as the Arnoldi iteration used in the ARPACK software [3]).

Computing Manifold Harmonics with Geogram

In Geogram, manifold harmonics are implemented in geogram/mesh/mesh_manifold_harmonics.h.

    void mesh_compute_manifold_harmonics(
	Mesh& M, index_t nb_eigens,
	LaplaceBeltramiDiscretization discretization,
	const std::string& attribute_name = "eigen",
	double shift = 0.0,
	bool print_spectrum = false
    );

The parameter discretization is one of:

  • COMBINATORIAL: 1.0 everywhere (graph Laplacian)
  • UNIFORM: combinatorial divided by node degree
  • FEM_P1: linear finite elements
  • FEM_P1_LUMPED: linear finite elements with lumped mass matrix (diagonal)

Under the hood, mesh_compute_manifold_harmonics() uses the shift_invert spectral transform, that makes the eigensolver dramatically faster: the eigensolver is faster for computing the eigenvectors associated with large eigenvalues, and unfortunately we want the other side of the spectrum. But it is easy to check that if a matrix M is non-singular, then M and M^{-1} have the same eigenvectors and the associated (non-zero) eigenvalues are the inverse. Note that one does not need to invert M explicitly (which would be prohibitive because the invert of a sparse matrix is not sparse in general). One only needs in practice to factor M instead. In addition, one can explore an arbitrary portion of the spectrum thanks to the parameter shift, that indicates from which part of the spectrum eigenvalues are searched.

When the number of eigenvalues/eigenfunctions is very large, not only the eigensolver may take a very long time to converge, but also storing the eigenvectors may be prohibitive. For this reason, there is another function, that computes the eigenvectors band by band (also by exploting the shift-invert spectral transform):

    void mesh_compute_manifold_harmonics_by_bands(
	Mesh& M, index_t nb_eigens,
	LaplaceBeltramiDiscretization discretization,
	ManifoldHarmonicsCallback callback,
	index_t nb_eigens_per_band = 30,
	double initial_shift = 0.0,
	void* client_data = nullptr
    );

For each eigenvalue/eigenvector, it will call the user-specified function callback that has the following prototype:

    typedef void (*ManifoldHarmonicsCallback)(
	index_t eigen_index,
	double eigen_val, const double* eigen_vector,
	void* client_data
    );

(and the client_data specified to mesh_compute_manifold_harmonics_by_bands() is passed back to it).

References

[1] Laplace-Beltrami Eigenfunctions Towards an Algorithm That "Understands" Geometry, Bruno Lévy, Shape Modeling International, 2006, HAL

[2] Spectral Geometry Processing with Manifold Harmonics, Bruno Vallet and Bruno Lévy, Computer Graphics Forum, 2008, HAL

[3] ARPACK software, homepage

Clone this wiki locally