/
init_surface.m
180 lines (147 loc) · 5.7 KB
/
init_surface.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
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
% init_surface Initializes a surface representation structure
%
% Usage:
%
% [surface] = init_surface (TRIV, X, Y, Z, [opt])
% [surface] = init_surface (surface, [opt])
%
% Description:
%
% Creates an internal representation of a surface for subsequent use
% by the TOSCA functions.
%
% Input:
%
% TRIV - ntx3 triangulation matrix with 1-based indices (as the one
% returned by the MATLAB function delaunay).
% X,Y,Z - vectors with nv vertex coordinates.
% surface - alternative way to specify the mesh as a struct, having .TRIV,
% .X, .Y, and .Z as its fields.
% opt - (optional) settings
% .verbose - Verbosity level (default: 1)
% 0 - no display
% 1 - display progress and surface parameters
%
% Output:
%
% surface - surface representation structure, comprising the following
% fields:
% .X,.Y,.Z - vertex coordinates
% .E - list of edges, each row corresponding to a pair of vertices
% defining an edge
% .TRIV - triangulation by vertex matrix, each row corresponds to a
% triangle represented by three vertices [v1 v2 v3]
% .TRIE - triangulation by edge matrix, each row corresponds to a
% triangle represented by three edges [e1 e2 e3]
% .ETRI - list of edge sharings, each row represents two triangle
% indices sharing the same edge from .E; one of the entries
% set to NaN denotes a boundary edge.
% .VTRI - cell array containing the list of triangles sharing a
% mesh vertex.
% .ADJV - cell array containing the list of adjacent vertices for
% each mesh vertex.
% .D - full matrix of geodesic distances between the points
% .diam - the mesh diameter (maximum geodesic distance)
% .nv - number of vertices
% .nt - number of triangles
% .ne - number of edges
% .genus - mesh genus
%
% References:
%
% [1] A. M. Bronstein, M. M. Bronstein, R. Kimmel, "Generalized multidimensional
% scaling: a framework for isometry-invariant partial surface matching",
% Proc. National Academy of Sciences (PNAS), Vol. 103/5, pp. 1168-1172,
% January 2006.
%
% [2] A. M. Bronstein, M. M. Bronstein, R. Kimmel, "Efficient computation of
% isometry-invariant distances between surfaces", SIAM Journal of
% Scientific Computing, Vol. 28/5, pp. 1812-1836, 2006.
%
% [3] A. M. Bronstein, M. M. Bronstein, R. Kimmel, "Calculus of non-rigid
% surfaces for geometry and texture manipulation", IEEE Trans.
% Visualization and Computer Graphics.
%
% TOSCA = Toolbox for Surface Comparison and Analysis
% Web: http://tosca.cs.technion.ac.il
% Version: 0.9
%
% (C) Copyright Alex Bronstein, 2007. All rights reserved.
%
% License:
%
% ANY ACADEMIC USE OF THIS CODE MUST CITE THE ABOVE REFERENCES.
% ANY COMMERCIAL USE PROHIBITED. PLEASE CONTACT THE AUTHORS FOR
% LICENSING TERMS. PROTECTED BY INTERNATIONAL INTELLECTUAL PROPERTY
% LAWS AND PATENTS PENDING.
function surface = init_surface (TRIV, X, Y, Z, opt)
% Constants
MAX_VERT_COUNT = 40000;
if nargin >= 4,
surface = [];
surface.X = X;
surface.Y = Y;
surface.Z = Z;
surface.TRIV = TRIV;
if nargin < 5, opt = []; end
else
surface = TRIV;
X = surface.X;
Y = surface.Y;
Z = surface.Z;
TRIV = surface.TRIV;
if nargin < 2, opt = []; else opt = X; end
end
% Options
if ~isfield(opt,'verbose'), opt.verbose = 1; else verbose = opt.verbose; end
if length(X) > MAX_VERT_COUNT,
error(sprintf('Current version of TOSCA can only work with meshes having less than %d vertices.\nUse remesh to reduce the vertex count.\n', MAX_VERT_COUNT));
end
% Initialize mesh structures
if opt.verbose > 0, fprintf(1,'Surface initialization\n'); end
if opt.verbose > 0, fprintf(1,'Computing mesh structures...\n'); end
% Reset mesh structures
surface.TRIE = [];
surface.ETRI = [];
% Create the edges structure
E = sort([TRIV(:,2) TRIV(:,3); TRIV(:,1) TRIV(:,3); TRIV(:,1) TRIV(:,2)],2);
[E,i,j] = unique(E,'rows');
surface.E = E;
% Create the triangle-by-edge structure
surface.TRIE = [j(1:end/3) j(end/3+1:2*end/3) j(2*end/3+1:end)];
% Create neighbouring triangles structure for each edge
% TO DO: more efficient C implementation!
surface.ETRI = zeros(length(E),2);
for e=1:length(E),
idx = find(surface.TRIE(:,1) == e | surface.TRIE(:,2) == e | surface.TRIE(:,3) == e)';
if length(idx(:))==1,
idx = [idx(:); NaN];
end
idx = idx(1:2);
surface.ETRI(e,:) = idx(:)';
end
% Create neighboring triangles structure for each vertex
% TO DO: more efficient C implementation!
surface.VTRI = cell(length(surface.X),1);
surface.ADJV = cell(length(surface.X),1);
for v = 1:length(surface.X),
surface.VTRI{v} = find(surface.TRIV(:,1) == v | surface.TRIV(:,2) == v | surface.TRIV(:,3) == v);
surface.ADJV{v} = setdiff(unique(surface.TRIV(surface.VTRI{v},:)), v);
end
% Initialize embedding structures
if opt.verbose > 0, fprintf(1,'Computing intrinsic geometry...\n'); end
% Compute geodesic distances
surface.D = fastmarch (surface);
surface.diam = max(surface.D(:));
% Compute mesh quantities
surface.nv = length(X);
surface.nt = size(TRIV,1);
surface.ne = length(E);
surface.genus = 1-0.5*(surface.nv-surface.ne+surface.nt);
if opt.verbose > 0,
fprintf(1,'Vertices\t\t%d\n', surface.nv);
fprintf(1,'Faces \t\t%d\n', surface.nt);
fprintf(1,'Edges \t\t%d\n', surface.ne);
fprintf(1,'Genus \t\t%.1f\n', surface.genus);
fprintf(1,'Diameter\t\t%d\n\n', surface.diam);
end