compute_greedy_homotopy_basis

Compute a greedy homotopy group basis based on the algorithm in paper[1]. Works for closed surface.

Two graph algorithms are needed: minimum spanning tree and shortest path, provided in Matlab's bioinformatics toolbox. We supply alternatives for these two functions, implemented purely in Matlab, which are a little slower and will be invoked when Maltab's built-in functions are not available.

  1. Erickson, Jeff, and Kim Whittlesey. "Greedy optimal homotopy and homology generators." Proceedings of the sixteenth annual ACM-SIAM symposium on Discrete algorithms. Society for Industrial and Applied Mathematics, 2005.

Contents

Syntax

hb = compute_greedy_homotopy_basis(face,vertex,bi)

Description

face  : double array, nf x 3, connectivity of mesh
vertex: double array, nv x 3, vertex of mesh
bi    : integer scaler, base point of homotopy group
hb: cell array, n x 1, a basis of homotopy group, each cell is a closed
    loop based at bi. Return empty for genus zero surface.

Contribution

Author : Wen Cheng Feng
Created: 2013/02/23
Revised: 2014/03/22 by Wen, optimize code with vectorized operation,
         remove dependency on Matlab's built-in function.
Revised: 2014/03/24 by Wen, add doc
Copyright 2014 Computational Geometry Group
Department of Mathematics, CUHK
http://www.math.cuhk.edu.hk/~lmlui
function hb = compute_greedy_homotopy_basis(face,vertex,bi)
% greedy_homotopy_basis(face,vertex,i) compute a homotopy
% basis of a high genus surface S = (face,vertex), at a basis point bi;
nf = size(face,1);
nv = size(vertex,1);
[edge,eif] = compute_edge(face);
[am,amd] = compute_adjacency_matrix(face);

% G is adjacency matrix constructed from the triangle mesh, with weight the
% edge length.
[I,J] = find(am);
el = sqrt(dot(vertex(I,:)-vertex(J,:),vertex(I,:)-vertex(J,:),2));
G = sparse(I,J,el,nv,nv);
G = (G + G'); % G is undirected

% the shortest path in graph G, with source node bi, T = G_pred is a the tree
if exist('graphshortestpath')
    [dist,path,pred] = graphshortestpath(G,bi,'METHOD','Dijkstra');
else
    [dist,path,pred] = dijkstra(G,bi);
end
% we always use column array
dist = dist(:);
pred = pred(:);
% amf is the dual graph of G
amf = compute_dual_graph(face);

% (G\T)*
% dual graph G* that do not consist edge correspond to edge in T = pred
I = (1:nv)';
I(bi) = [];
J = pred(I);
% (I2,J2) and (J2,I2) are faces corresponding to edge (I,pred(I))
I2 = full(amd(I+(J-1)*nv));
J2 = full(amd(J+(I-1)*nv));
ind = (I2==0 | J2==0);
I2(ind) = [];
J2(ind) = [];
amf(I2+(J2-1)*nf) = 0;
amf(J2+(I2-1)*nf) = 0;

% tree is the maximum spanning tree of (G\T)*, where the weight of any
% edge e* is length(shortest_loop(e))
[I,J] = find(amf);
ind = (eif(:,1)==-1 | eif(:,2)==-1);
eif(ind,:) = [];
edge(ind,:) = [];
F2E = sparse([eif(:,1);eif(:,2)],[eif(:,2);eif(:,1)],[edge(:,1);edge(:,2)],nf,nf);
ei = [F2E(I+(J-1)*nf),F2E(J+(I-1)*nf)];
dvi = vertex(ei(:,1),:)-vertex(ei(:,2),:);
V = -(dist(ei(:,1))+dist(ei(:,2))+sqrt(dot(dvi,dvi,2)));
amf_w = sparse(I,J,V,nf,nf);
if exist('graphminspantree')
    tree = graphminspantree(amf_w,'METHOD','Kruskal');
else
    tree = minimum_spanning_tree(amf_w);
end

% G2 is the graph, with edges neither in T nor are crossed by edges in T*
G2 = G;
I = (1:nv)';
I(bi) = [];
J = pred(I);
% (I,J) are edges in T
G2(I+(J-1)*nv) = 0;
G2(J+(I-1)*nv) = 0;

[I,J] = find(tree);
% ei are edges in T*
ei = [F2E(I+(J-1)*nf),F2E(J+(I-1)*nf)];
index = [ei(:,1)+(ei(:,2)-1)*nv;ei(:,2)+(ei(:,1)-1)*nv];
G2(index) = 0;

% the greedy homotopy basis consists of all loops (e), where e is an edge
%  of G2.
G2 = tril(G2);
[I,J] = find(G2);
hb = cell(size(I));
for i = 1:length(I)
    pi = path{I(i)};
    pj = path{J(i)};
    hb{i} = [pi(:);flipud(pj(:))];
end
if isempty(I)
    hb = [];
end