Compute principal curvatures on mesh vertices.
[C1, C2] = meshCurvatures(VERTICES, FACES)
Computes the principal curvatures C1 and C2 for each vertex of the mesh
defined by VERTICES and FACES.
[C1, C2] = meshCurvatures(..., PNAME, PVALUE)
Provides additional input arguments based on a list of name-value pairs
of arguments. Parameter names can be:
* 'SmoothingSteps' (integer, default: 3)
Specifies the number of steps for smoothing vertex curvature
tensors.
* 'Verbose' (boolean, default: true)
Displays details about algorithm processing.
* 'ShowProgress' (boolean, default: true)
Displays a text-based progress bar.
Algorithm
The function is adapted from the "compute_curvature" function, in the
"toolbox_graph" fro Gabriel Peyre.
The basic idea is to define a curvature tensor for each edge, by
assigning a minimum curvature equal to zero along the edge, and a
maximum curvature equal to the dihedral angle across the edge.
Averaging around the neighbors of a vertex v yields a summation formula
over the neighbor edges to compute the curvature tensor of a vertex:
1
C(v) = ---- Sum \beta(e) || e \cap A(v) || ebar ebar^t
A(v) {e \in A(v)}
where:
* A(v) is the neighborhood region, usually defined as a 'ring' around
the vertex v
* beta(e) is the dihedral angle between the normals of the two faces
incident to edge e
* || e \cap A(v) || is the length of e (more exactly, the length of the
part of e contained within the neighborhood region
* ebar is the normalized edge
The curvature tensor is then decomposed into C = P D P^-1, with P
containing main direction vectors and normal, and D being a diagonal
matrix with the two main curvatures and zero along the diagonal.
References
* David Cohen-Steiner and Jean-Marie Morvan (2003).
"Restricted Delaunay triangulations and normal cycle".
In Proc. 19th Annual ACM Symposium on Computational Geometry,
pages 237-246.
* Pierre Alliez, David Cohen-Steiner, Olivier Devillers, Bruno Levy,
and Mathieu Desbrun (2003). "Anisotropic Polygonal Remeshing".
ACM Transactions on Graphics.
(SIGGRAPH '2003 Conference Proceedings)
* Mario Botsch, Leif Kobbelt, M. Pauly, P. Alliez, B. Levy (2010).
"Polygon Mesh Processing", Taylor and Francis Group, New York.
Example
[v, f] = torusMesh;
f2 = triangulateFaces(f);
[c1, c2] = meshCurvatures(v, f2);
figure; hold on; axis equal; view(3);
drawMesh(v, f2, 'VertexColor', c1 .* c2);
See also
meshes3d, drawMesh, triangulateFaces0001 function [C1, C2, U1, U2, H, K, N] = meshCurvatures(vertices, faces, varargin) 0002 % Compute principal curvatures on mesh vertices. 0003 % 0004 % [C1, C2] = meshCurvatures(VERTICES, FACES) 0005 % Computes the principal curvatures C1 and C2 for each vertex of the mesh 0006 % defined by VERTICES and FACES. 0007 % 0008 % [C1, C2] = meshCurvatures(..., PNAME, PVALUE) 0009 % Provides additional input arguments based on a list of name-value pairs 0010 % of arguments. Parameter names can be: 0011 % * 'SmoothingSteps' (integer, default: 3) 0012 % Specifies the number of steps for smoothing vertex curvature 0013 % tensors. 0014 % * 'Verbose' (boolean, default: true) 0015 % Displays details about algorithm processing. 0016 % * 'ShowProgress' (boolean, default: true) 0017 % Displays a text-based progress bar. 0018 % 0019 % Algorithm 0020 % The function is adapted from the "compute_curvature" function, in the 0021 % "toolbox_graph" fro Gabriel Peyre. 0022 % The basic idea is to define a curvature tensor for each edge, by 0023 % assigning a minimum curvature equal to zero along the edge, and a 0024 % maximum curvature equal to the dihedral angle across the edge. 0025 % Averaging around the neighbors of a vertex v yields a summation formula 0026 % over the neighbor edges to compute the curvature tensor of a vertex: 0027 % 1 0028 % C(v) = ---- Sum \beta(e) || e \cap A(v) || ebar ebar^t 0029 % A(v) {e \in A(v)} 0030 % where: 0031 % * A(v) is the neighborhood region, usually defined as a 'ring' around 0032 % the vertex v 0033 % * beta(e) is the dihedral angle between the normals of the two faces 0034 % incident to edge e 0035 % * || e \cap A(v) || is the length of e (more exactly, the length of the 0036 % part of e contained within the neighborhood region 0037 % * ebar is the normalized edge 0038 % 0039 % The curvature tensor is then decomposed into C = P D P^-1, with P 0040 % containing main direction vectors and normal, and D being a diagonal 0041 % matrix with the two main curvatures and zero along the diagonal. 0042 % 0043 % References 0044 % * David Cohen-Steiner and Jean-Marie Morvan (2003). 0045 % "Restricted Delaunay triangulations and normal cycle". 0046 % In Proc. 19th Annual ACM Symposium on Computational Geometry, 0047 % pages 237-246. 0048 % * Pierre Alliez, David Cohen-Steiner, Olivier Devillers, Bruno Levy, 0049 % and Mathieu Desbrun (2003). "Anisotropic Polygonal Remeshing". 0050 % ACM Transactions on Graphics. 0051 % (SIGGRAPH '2003 Conference Proceedings) 0052 % * Mario Botsch, Leif Kobbelt, M. Pauly, P. Alliez, B. Levy (2010). 0053 % "Polygon Mesh Processing", Taylor and Francis Group, New York. 0054 % 0055 % Example 0056 % [v, f] = torusMesh; 0057 % f2 = triangulateFaces(f); 0058 % [c1, c2] = meshCurvatures(v, f2); 0059 % figure; hold on; axis equal; view(3); 0060 % drawMesh(v, f2, 'VertexColor', c1 .* c2); 0061 % 0062 % See also 0063 % meshes3d, drawMesh, triangulateFaces 0064 % 0065 0066 % ------ 0067 % Author: David Legland 0068 % e-mail: david.legland@inrae.fr 0069 % INRAE - BIA Research Unit - BIBS Platform (Nantes) 0070 % Created: 2021-09-21, using Matlab 9.10.0.1684407 (R2021a) Update 3 0071 % Copyright 2021 INRAE. 0072 0073 0074 %% Process input arguments 0075 0076 % default values for options 0077 nIters = 3; 0078 verbose = true; 0079 showProgress = true; 0080 0081 while length(varargin) > 1 0082 name = varargin{1}; 0083 if strcmpi(name, 'SmoothingSteps') 0084 nIters = varargin{2}; 0085 elseif strcmpi(name, 'Verbose') 0086 verbose = varargin{2}; 0087 elseif strcmpi(name, 'ShowProgress') 0088 showProgress = varargin{2}; 0089 else 0090 error('Unknown option: %s', name); 0091 end 0092 varargin(1:2) = []; 0093 end 0094 0095 % validate vertices 0096 if ~isnumeric(vertices) || size(vertices, 2) ~= 3 0097 error('Requires vertices to be a N-by-3 numeric array'); 0098 end 0099 0100 % ensure faces are triangular 0101 if ~isnumeric(faces) || size(faces, 2) > 3 0102 warning('requires triangle mesh, forces triangulation'); 0103 faces = triangulateFaces(faces); 0104 end 0105 0106 0107 %% Retrieve adjacency relationships 0108 0109 if verbose 0110 disp('compute adjacencies'); 0111 end 0112 0113 % number of elements of each type 0114 nv = size(vertices, 1); 0115 nf = size(faces, 1); 0116 0117 % ev1 and ev2 are indices of source and target vertex of each edge 0118 % (recomputed later) 0119 ev1 = [faces(:,1); faces(:,2); faces(:,3)]; 0120 ev2 = [faces(:,2); faces(:,3); faces(:,1)]; 0121 0122 % Compute sparse matrix representing edge-to-face adjacency 0123 s = [1:nf 1:nf 1:nf]'; 0124 A = sparse(ev1, ev2, s, nv, nv); 0125 0126 % converts sparse matrix to indices of adjacent vertices and faces 0127 [~, ~, ef1] = find(A); % index of 'right' face 0128 [ev1, ev2, ef2] = find(A'); % index of 'left' face, and of vertices 0129 0130 % edges are consdered twice (one for each vertex) 0131 % keep only the edge with lower source index 0132 inds = find(ev1 < ev2); 0133 ef1 = ef1(inds); 0134 ef2 = ef2(inds); 0135 ev1 = ev1(inds); 0136 ev2 = ev2(inds); 0137 0138 % number of edges 0139 ne = length(ev1); 0140 0141 0142 %% Compute geometry features 0143 0144 % compute edge direction vectors 0145 edgeVectors = vertices(ev2,:) - vertices(ev1,:); 0146 0147 % normalize edge direction vecotrs 0148 d = sqrt(sum(edgeVectors.^2, 2)); 0149 edgeVectors = bsxfun(@rdivide, edgeVectors, d); 0150 0151 % avoid too large numerics 0152 d = d ./ mean(d); 0153 0154 % normals to faces 0155 normals = meshFaceNormals(vertices, faces); 0156 0157 % ensure normals point outward the mesh 0158 if meshVolume(vertices, faces) < 0 0159 normals = -normals; 0160 end 0161 0162 % inner product of normals 0163 dp = sum(normals(ef1, :) .* normals(ef2, :), 2); 0164 0165 % compute the (unsigned) dihedral angle between the normals of the two 0166 % faces incident to each edge 0167 beta = acos(min(max(dp, -1), 1)); 0168 0169 % relatice orientation of face normals cross product and edge orientation 0170 cp = crossProduct3d(normals(ef1, :), normals(ef2, :)); 0171 si = sign(sum(cp .* edgeVectors, 2)); 0172 0173 % compute signed dihedral angle 0174 beta = beta .* si; 0175 0176 0177 %% Compute tensors 0178 0179 if verbose 0180 disp('compute edge tensors'); 0181 end 0182 0183 % curvature tensor of each edge 0184 T = zeros(3, 3, ne); 0185 for i = 1:3 0186 for j = 1:i 0187 T(i, j, :) = reshape(edgeVectors(:,i) .* edgeVectors(:,j), 1, 1, ne); 0188 T(j, i, :) = T(i, j, :); 0189 end 0190 end 0191 T = bsxfun(@times, T, reshape(d .* beta, [1 1 ne])); 0192 0193 % curvature tensor of each vertex by pooling edge tensors 0194 Tv = zeros(3, 3, nv); 0195 w = zeros(1, 1, nv); 0196 for k = 1:ne 0197 if showProgress 0198 displayProgress(k, ne); 0199 end 0200 Tv(:,:,ev1(k)) = Tv(:,:,ev1(k)) + T(:,:,k); 0201 Tv(:,:,ev2(k)) = Tv(:,:,ev2(k)) + T(:,:,k); 0202 w(:,:,ev1(k)) = w(:,:,ev1(k)) + 1; 0203 w(:,:,ev2(k)) = w(:,:,ev2(k)) + 1; 0204 end 0205 w(w < eps) = 1; 0206 Tv = Tv ./ repmat(w, [3 3 1]); 0207 0208 if verbose 0209 disp('average vertex tensors'); 0210 end 0211 0212 % apply smoothing on the tensor field 0213 for i = 1:3 0214 for j = 1:3 0215 a = Tv(i, j, :); 0216 a = smoothMeshFunction(vertices, faces, a(:), nIters); 0217 Tv(i, j, :) = reshape(a, [1 1 nv]); 0218 end 0219 end 0220 0221 0222 %% Retrieve curvatures and eigen vectors from tensors 0223 0224 if verbose 0225 disp('retrieve curvatures'); 0226 end 0227 0228 % allocate memory 0229 U = zeros(3, 3, nv); 0230 D = zeros(3, nv); 0231 0232 % iterate over vertices 0233 for k = 1:nv 0234 % display progress 0235 if showProgress 0236 displayProgress(k,nv); 0237 end 0238 0239 % extract eigenvectors and eigenvalues for current vertex 0240 [u, d] = eig(Tv(:,:,k)); 0241 d = real(diag(d)); 0242 0243 % sort acording to [normal, min curv, max curv] 0244 [~, I] = sort(abs(d)); 0245 D(:, k) = d(I); 0246 U(:, :, k) = real(u(:,I)); 0247 end 0248 0249 % retrieve main curvatures and associated directions 0250 C1 = D(2,:)'; 0251 C2 = D(3,:)'; 0252 U1 = squeeze(U(:,3,:))'; 0253 U2 = squeeze(U(:,2,:))'; 0254 0255 % enforce C1 < C2 0256 inds = find(C1 > C2); 0257 C1tmp = C1; 0258 U1tmp = U1; 0259 C1(inds) = C2(inds); 0260 C2(inds) = C1tmp(inds); 0261 U1(inds,:) = U2(inds,:); 0262 U2(inds,:) = U1tmp(inds,:); 0263 0264 % compute optional output arguments 0265 if nargout > 4 0266 % average and gaussian curvatures 0267 H = (C1 + C2) / 2; 0268 K = C1 .* C2; 0269 0270 if nargout > 6 0271 % normal vector for each vertex 0272 N = squeeze(U(:,1,:))'; 0273 end 0274 end 0275 0276 0277 function displayProgress(n, N) 0278 % Display the progress of current step using a text-based progress bar. 0279 % 0280 % based on the 'progressbar' function in G. Peyre's Graph Toolbox. 0281 0282 % width of the progress bar 0283 w = 20; 0284 0285 % compute progress ratio as an integer betsween 0 and w 0286 p = min( floor(n/N*(w+1)), w); 0287 0288 global pprev; 0289 if isempty(pprev) 0290 pprev = -1; 0291 end 0292 0293 if p ~= pprev 0294 str1 = repmat('*', 1, p); 0295 str2 = repmat('.', 1, w-p); 0296 str = sprintf('[%s%s]', str1, str2); 0297 if n > 1 0298 % clear previous string 0299 fprintf(repmat('\b', [1 length(str)])); 0300 end 0301 fprintf(str); 0302 end 0303 0304 pprev = p; 0305 if n == N 0306 fprintf('\n'); 0307 end