# Toolbox graph - A toolbox to process graph and triangulated meshes

Copyright (c) 2008 Gabriel Peyre

## Contents

- 3D Mesh Loading and Displaying
- Differential Calculus on 3D Mesh
- Mesh Filtering
- Mesh parameterization
- Tetrahedral Mesh Loading and Displaying

The toolbox can be downloaded from Matlab Central http://www.mathworks.com/matlabcentral/fileexchange/loadFile.do?objectId=5355&objectType=FILE

First add to the path additional scripts.

    path(path, 'toolbox/');
    path(path, 'off/');
    clear options;

## 3D Mesh Loading and Displaying

A 3D mesh is composed of a `vertex` array of size `(3,n)` that contains the position in 3D (or sometimes 2D) of the vertices, and a `face` array of dimension `(3,m)` that contains the indexes of each triangulated face.

    % load the mesh
    name = 'elephant-50kv';
    options.name = name; % useful for displaying
    [vertex,faces] = read_mesh(name);
    % display the mesh
    clf;
    plot_mesh(vertex, faces);
    shading interp;

You can zoom on the mesh and display its triangulated faces

    clf;
    for i=1:4
        subplot(2,2,i);
        plot_mesh(vertex, faces);
        shading faceted;
        zoom(1.8^(i+1));
    end

You can compute the normal to the mesh, and displace the position of the vertices along the normal.

    options.face_vertex_color =  [];
    vertex1 = perform_normal_displacement(vertex,faces,.03);

    clf;
    subplot(1,2,1);
    plot_mesh(vertex,faces,options); shading interp; axis tight;
    subplot(1,2,2);
    plot_mesh(vertex1,faces,options); shading interp; axis tight;

## Differential Calculus on 3D Mesh

You can compute matrix (sparse) of differential and averaging operators on 3D meshes. You can also compute curvature and normals.

You can compute normals and display them.

    % load a mesh
    name = 'mushroom';
    options.name = name; % useful for displaying
    [vertex,faces] = read_mesh(name);
    % compute normal per vertex and per face
    [normal,normalf] = compute_normal(vertex,faces);
    % display
    options.normal = normal;
    clf; plot_mesh(vertex,faces,options); shading interp; axis tight;
    options.normal = [];

You can compute curvature quantities using local avergaging and PCA analysis.

    % load the mesh
    name = 'elephant-50kv';
    options.name = name; % useful for displaying
    [vertex,faces] = read_mesh(name);
    % compute the curvature
    options.curvature_smoothing = 10;
    options.verb = 0;
    [Umin,Umax,Cmin,Cmax,Cmean,Cgauss,Normal] = compute_curvature(vertex,faces,options);
    % display
    clf;
    subplot(1,2,1);
    options.face_vertex_color = perform_saturation(Cgauss,1.2);
    plot_mesh(vertex,faces, options); shading interp; colormap jet(256);
    title('Gaussian curvature');
    subplot(1,2,2);
    options.face_vertex_color = perform_saturation(abs(Cmin)+abs(Cmax),1.2);
    plot_mesh(vertex,faces, options); shading interp; colormap jet(256);
    title('Total curvature');

A laplacian is an high pass operator that compute second order derivatives.

    name = 'elephant-50kv';
    options.name = name; % useful for displaying
    [vertex,faces] = read_mesh(name);
    % options for the Laplacians
    laplacian_type = 'conformal'; % slow to compute
    laplacian_type = 'combinatorial'; % fast, but inexact
    laplacian_type = 'distance'; % fast and more accurate
    % compute a normalized Laplacian
    options.symmetrize = 0;
    options.normalize = 1;
    L = compute_mesh_laplacian(vertex,faces,laplacian_type,options);
    % compute a non normalized Laplacian
    options.symmetrize = 0;
    options.normalize = 1;
    L0 = compute_mesh_laplacian(vertex,faces,laplacian_type,options);

Compute and display some eigenvectors of the Laplacian. They are the equivalent of the Fourier sinusoids on surfaces.

    % compute a few eigenmodes
    nb = 60;
    opts.disp = 0;
    [V,D] = eigs(L0,nb,'SM',opts); % warning : it takes lot of time
    V = real(V(:,end:-1:1));
    % display them on the mesh
    ilist = round(linspace(3,nb-3, 6));
    tau=2.2; % saturation for display
    clf;
    for i=1:length(ilist)
        % subplot(1,length(ilist),i);
        v = real(V(:,ilist(i)));
        v = clamp( v/std(v),-tau,tau );
        options.face_vertex_color = v;
        subplot(2,3,i);
        plot_mesh(vertex,faces,options);
        shading interp; camlight; axis tight; % zoom(zoomf);
        colormap jet(256);
    end

## Mesh Filtering

You can filter a mesh by applying a filter operator, or iterating several time the filtering, which corresponds to the resolution of the heat equation on the mesh

First we create a noisy mesh by displacement of the vertices along the normal direction (those are the most distructive displacements).

    normals = compute_normal(vertex,faces);
    rho = randn(1,size(vertex,2))*.02;
    vertex1 = vertex + repmat(rho,3,1).*normals;

Then we compute a filtering operator (low pass)

    % parameters for the operator
    laplacian_type = 'distance';
    options.symmetrize = 0;
    options.normalize = 1; % it must be normalized for filtering
    options.verb = 0;
    W = compute_mesh_weight(vertex,faces,laplacian_type,options);
    % This is the corresponding laplacian
    L = compute_mesh_laplacian(vertex,faces,laplacian_type,options);

We apply iteratively the operator to get filtered versions of the mesh.

    vertex2 = vertex1;
    clf;
    options.face_vertex_color = [];
    for i=1:6
        subplot(2,3,i);
        plot_mesh(vertex2,faces,options); axis tight; shading interp;
        vertex2 = (W*(W*vertex2'))';
    end

Something better is to update `W` at each filtering step. You can also try to enforce various quantity during the iterations, for instance the volume of the mesh.

Instead of iterative filtering, one can solve for the heat diffusion flow.

    Tlist = [0 2 5 10 20 40]/4;
    options.dt = 0.3;
    options.face_vertex_color = [];
    clf;
    for k=1:length(Tlist)
        options.Tmax = Tlist(k);
        vertex2 = perform_mesh_heat_diffusion(vertex1,faces,L,options);
        subplot(2,3,k);
        plot_mesh(vertex2,faces,options);
        shading interp; axis tight;
    end

You can try to select the optimal number of iterations for applying `W` or the optimal time for the heat equation so that to minimize the denoising error.

## Mesh parameterization

Mesh parameterization corresponds to the computation of a 2D position for each vertex of a 3D mesh, so that the 2D->3D mapping is bijective.

First load a small mesh and set up parameters

    name = 'nefertiti'; options.name = name;
    [vertex,faces] = read_mesh(name);
    A = triangulation2adjacency(faces);
    % you can try with other boundary type
    options.boundary = 'circle';
    % you can try with other Laplacians
    options.laplacian = 'conformal';

The fixed boundary parameterization is guaranteed to be valid is the mesh is homeomorphic to a disk and if the boundary of the planar domain is convex. This is the Tutte embedding theorem.

    % compute the layout in 2D
    options.method = 'parameterization';
    options.verb = 0;
    xy = compute_parameterization(vertex,faces,options);
    % display the parameterization
    clf;
    subplot(1,2,1);
    plot_mesh(vertex,faces,options); shading faceted; axis tight;
    subplot(1,2,2);
    plot_graph(A,xy,'k.-'); axis tight;
    title('Parameterization');

In order to reduct the distortion, it is possible to relax the boundary condition and sovle for a Neumann boundary condition. The parameterization is not guaranteed to be valid anymore.

    % compute the parameterization
    options.method = 'freeboundary';
    xy1 = compute_parameterization(vertex,faces,options);
    % display the parameterization
    clf;
    subplot(1,2,1);
    plot_graph(A,xy,'k.-'); axis tight;
    title('Fixed boundary');
    subplot(1,2,2);
    plot_graph(A,xy1,'k.-'); axis tight;
    title('Free boundary');

Another option is to use dimensionality reduction technics, for instance Laplacian eigenmaps (first eigenvectors of the laplacian) or Isomap (geodesic distances computations).

    % perform dimensionality reduction
    options.method = 'flattening';
    xy2 = compute_parameterization(vertex,faces,options);
    options.method = 'isomap';
    xy3 = compute_parameterization(vertex,faces,options);
    % display
    clf;
    subplot(1,2,1);
    plot_graph(A,xy2,'k.-'); axis tight;
    title('Laplacian eigenmaps');
    subplot(1,2,2);
    plot_graph(A,xy3,'k.-'); axis tight;
    title('Isomap');

## Tetrahedral Mesh Loading and Displaying

You can load and display volumetric tetrahedral meshes. Important: .tet files and not included in the toolbox distribution (too large files). You should download them from

http://www.aimatshape.net/

    [vertex,face] = read_tet('../toolbox_graph_data/tet/hand.tet');
    clf;
    plot_mesh(vertex,face,options);