laplace_beltrami

Laplace Beltrami operator on the mesh.

Cotangent formula is used, while there are some variants:

For comparison and convergence analysis, see paper [4]

  1. K. Polthier. Computational Aspects of Discrete Minimal Surfaces. In Proc. of the Clay Summer School on Global Theory of Minimal Surfaces, J. Hass, D. Hoffman, A. Jaffe, H. Rosenberg, R. Schoen, M. Wolf (Eds.), to appear, 2002.
  2. M. Meyer, M. Desbrun, P. Schröder, and A. Barr. Discrete Differential-Geometry Operator for Triangulated 2-manifolds. In Proc. VisMath'02, Berlin, Germany, 2002.
  3. M. Desbrun, M. Meyer, P. Schröder, and A. H. Barr. Implicit Fairing of Irregular Meshes using Diffusion and Curvature Flow. SIGGRAPH99, pages 317-324, 1999.
  4. Xu, Guoliang. "Convergent discrete laplace-beltrami operators over triangular surfaces." Geometric Modeling and Processing, 2004. Proceedings. IEEE, 2004.

Contents

Syntax

A = laplace_beltrami(face,vertex)
A = laplace_beltrami(face,vertex,method)

Description

face  : double array, nf x 3, connectivity of mesh
vertex: double array, nv x 3, vertex of mesh
method: string, optional, method of cotangent formula, can be one of
        three: 'Polthier', 'Meyer', 'Desbrun'. Default is 'Polthier'.
A: sparse matrix, nv x nv, Laplace Beltrami operator

Example

A = laplace_beltrami(face,vertex)
A = laplace_beltrami(face,vertex,'Polthier') % same as last
A = laplace_beltrami(face,vertex,'Meyer')
A = laplace_beltrami(face,vertex,'Desbrun')

Contribution

Author : Wen Cheng Feng
Created: 2014/03/03
Revised: 2014/03/03 by Wen, add more cotangent formula variants, not
         implemented
Revised: 2014/03/23 by Wen, add doc
Revised: 2014/03/28 by Wen, add code for 'Meyer' and 'Desbrun' methods
Copyright 2014 Computational Geometry Group
Department of Mathematics, CUHK
http://www.math.cuhk.edu.hk/~lmlui
function A = laplace_beltrami(face,vertex,method)
% default method is 'Polthier'
if nargin == 2
    method = 'Polthier';
end
[edge,eif] = compute_edge(face);
ne = size(edge,1);
ew = zeros(ne,1);
ind = eif(:,1)>0;
ev1 = sum(face(eif(ind,1),:),2) - sum(edge(ind,:),2);
ct1 = cot2(vertex(ev1,:),vertex(edge(ind,1),:),vertex(edge(ind,2),:));
ew(ind) = ew(ind) + ct1;
ind = eif(:,2)>0;
ev2 = sum(face(eif(ind,2),:),2) - sum(edge(ind,:),2);
ct2 = cot2(vertex(ev2,:),vertex(edge(ind,1),:),vertex(edge(ind,2),:));
ew(ind) = ew(ind) + ct2;

switch method
    case 'Polthier';
        A = sparse([edge(:,1);edge(:,2)],[edge(:,2);edge(:,1)],[ew;ew]/2);
        sA = sum(A,2);
        A = A - diag(sA);
    case 'Meyer'
        va = vertex_area(face,vertex,'mixed');
        ew = (ew./va(edge(:,1))+ew./va(edge(:,2)))/2;
        A = sparse([edge(:,1);edge(:,2)],[edge(:,2);edge(:,1)],[ew;ew]);
        sA = sum(A,2);
        A = A - diag(sA);
    case 'Desbrun'
        va = vertex_area(face,vertex,'one_ring');
        ew = (ew./va(edge(:,1))+ew./va(edge(:,2)))/2*3;
        A = sparse([edge(:,1);edge(:,2)],[edge(:,2);edge(:,1)],[ew;ew]);
        sA = sum(A,2);
        A = A - diag(sA);
    otherwise
        error('Wrong method. Available methods are: Polthier,Meyer,Desbrun.')
end

function ct = cot2(pi,pj,pk)
a = sqrt(dot(pj-pk,pj-pk,2));
b = sqrt(dot(pk-pi,pk-pi,2));
c = sqrt(dot(pi-pj,pi-pj,2));
cs = (b.*b+c.*c-a.*a)./(2.*b.*c);
ss2 = 1-cs.*cs;
ss2(ss2<0) = 0;
ss2(ss2>1) = 1;
ss = sqrt(ss2);
ct = cs./ss;