-
Notifications
You must be signed in to change notification settings - Fork 5
/
Copy pathmain_qslim_anisotropic.m
113 lines (97 loc) · 3.05 KB
/
main_qslim_anisotropic.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
clear all; close all;
addpath(genpath('./gptoolbox'))
addpath(genpath('./utils'))
addpath(genpath('./comparison'))
% coarsening parameter
numVc = 600;
numEig = numVc-1;
numEig_vis = 100; % number of eigs to visualize
normalizationScale = sqrt(numVc);
numRings = 1; % number of rings (change this to set 1-, 2- and 3-ring)
%% read mesh
objname = './data/models/armadillo.obj';
% objname = './data/models/spot.obj';
[filepath,name,ext] = fileparts(objname);
if strcmp(ext, '.obj')
[V,F] = readOBJ(objname);
elseif strcmp(ext, '.stl')
[V,F] = readSTL(objname);
elseif strcmp(ext, '.ply')
[V,F] = readPLY(objname);
end
V = normalizationScale * normalizeUnitArea(V,F);
%% generate coarse mesh using [Garland & Heckbert 1997]
[Vc,Fc,R] = QSlim_withR(V,F,numVc);
fprintf('============Finish generating coarse mesh connection using [Garland & Heckbert 1997]==========\n');
fprintf('==================================Our algorithm starts========================================\n');
%% compute cotangent laplacian and mass matrix
alpha = 50.0;
options.n_eigenvalues = numEig;
L = -anisoLaplace(V,F,alpha,options);
M = massmatrix(V,F);
Lc = -anisoLaplace(Vc,Fc,alpha,options);
Mc = massmatrix(Vc,Fc);
% set the sparsity pattern
Xsp = (Lc~=0);
Xsp = (Xsp^numRings ~= 0); % sparsity pattern
X_init = Lc;
%% ADMM solve with chordal decomposition
[X, state] = chordal_coarsening(L, M, R, Xsp, Mc, X_init, 'numEig', numEig);
fprintf("our solve time: %.4e\n", state.time_admm);
%% visualize the results
% plot results
[eVal,eVec] = eigsReal(L, M, numEig_vis);
[eVal_ours, eVec_ours] = eigsReal(X, Mc, numEig_vis);
[eVal_cotan, eVec_cotan] = eigsReal(Lc, Mc, numEig_vis);
% make sure the eigenvectors are mass-orthogonal
eVec_ours = massOrthogonal(eVec_ours, Mc);
eVec = massOrthogonal(eVec, M);
eVec_cotan = massOrthogonal(eVec_cotan, Mc);
fMap_ours = eVec_ours' * Mc * R * eVec;
fMap_cotan = eVec_cotan' * Mc * R * eVec;
figure(1)
subplot(1,2,1)
plotFMap(fMap_ours(1:numEig_vis, 1:numEig_vis))
title('fmap (ours)')
subplot(1,2,2)
plotFMap(fMap_cotan(1:numEig_vis, 1:numEig_vis))
title('fmap cotan Laplace')
figure(2)
plot(log10(state.objHis))
xlabel('iterations')
ylabel('log10 energy value')
figure(3)
plot(eVal)
hold on
plot(eVal_ours)
plot(eVal_cotan)
legend('original', 'ours','cotan Laplace')
ylabel('eigenvalues')
% visualize one eigenfunctions (this may have sign flip)
figure(4)
subplot(1,3,1)
plotMesh(V,F,eVec(:,20))
title('original');
subplot(1,3,2)
plotMesh(Vc,Fc,eVec_ours(:,20))
title('ours');
subplot(1,3,3)
plotMesh(Vc,Fc,eVec_cotan(:,20))
title('[Garland & Heckbert 1997]');
axis equal off
sgtitle('visualize one eigen function (may have sign flip)')
figure(5)
subplot(1,2,1)
spy(Lc);
title('[Garland & Heckbert 1997]');
subplot(1,2,2)
spy(X)
title('our pattern');
% visualize functional map
fmap = fmap_operator(numEig_vis, L, M);
other = struct();
[other.C, other.nC, other.nL, other.nD] = fmap(R, Lc, Mc);
fig_fmap(6, '[Garland & Heckbert 1997]', other.C, other.nL, other.nD);
our = struct();
[our.C, our.nC, our.nL, our.nD] = fmap(R, X, Mc);
fig_fmap(7, '[Ours]', our.C, our.nL, our.nD);