MESHCURVATURES 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" from 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 %MESHCURVATURES 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" from 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 % Created: 2021-09-21, using Matlab 9.10.0.1684407 (R2021a) Update 3 0070 % Copyright 2021-2024 INRAE - BIA Research Unit - BIBS Platform (Nantes) 0071 0072 %% Process input arguments 0073 0074 % default values for options 0075 nIters = 3; 0076 verbose = true; 0077 showProgress = true; 0078 0079 while length(varargin) > 1 0080 name = varargin{1}; 0081 if strcmpi(name, 'SmoothingSteps') 0082 nIters = varargin{2}; 0083 elseif strcmpi(name, 'Verbose') 0084 verbose = varargin{2}; 0085 elseif strcmpi(name, 'ShowProgress') 0086 showProgress = varargin{2}; 0087 else 0088 error('Unknown option: %s', name); 0089 end 0090 varargin(1:2) = []; 0091 end 0092 0093 % validate vertices 0094 if ~isnumeric(vertices) || size(vertices, 2) ~= 3 0095 error('Requires vertices to be a N-by-3 numeric array'); 0096 end 0097 0098 % ensure faces are triangular 0099 if ~isnumeric(faces) || size(faces, 2) > 3 0100 warning('requires triangle mesh, forces triangulation'); 0101 faces = triangulateFaces(faces); 0102 end 0103 0104 0105 %% Retrieve adjacency relationships 0106 0107 if verbose 0108 disp('compute adjacencies'); 0109 end 0110 0111 % number of elements of each type 0112 nv = size(vertices, 1); 0113 nf = size(faces, 1); 0114 0115 % ev1 and ev2 are indices of source and target vertex of each edge 0116 % (recomputed later) 0117 ev1 = [faces(:,1); faces(:,2); faces(:,3)]; 0118 ev2 = [faces(:,2); faces(:,3); faces(:,1)]; 0119 0120 % Compute sparse matrix representing edge-to-face adjacency 0121 s = [1:nf 1:nf 1:nf]'; 0122 A = sparse(ev1, ev2, s, nv, nv); 0123 0124 % converts sparse matrix to indices of adjacent vertices and faces 0125 [~, ~, ef1] = find(A); % index of 'right' face 0126 [ev1, ev2, ef2] = find(A'); % index of 'left' face, and of vertices 0127 0128 % edges are consdered twice (one for each vertex) 0129 % keep only the edge with lower source index 0130 inds = find(ev1 < ev2); 0131 ef1 = ef1(inds); 0132 ef2 = ef2(inds); 0133 ev1 = ev1(inds); 0134 ev2 = ev2(inds); 0135 0136 % number of edges 0137 ne = length(ev1); 0138 0139 0140 %% Compute geometry features 0141 0142 % compute edge direction vectors 0143 edgeVectors = vertices(ev2,:) - vertices(ev1,:); 0144 0145 % normalize edge direction vecotrs 0146 d = sqrt(sum(edgeVectors.^2, 2)); 0147 edgeVectors = bsxfun(@rdivide, edgeVectors, d); 0148 0149 % avoid too large numerics 0150 d = d ./ mean(d); 0151 0152 % normals to faces 0153 normals = meshFaceNormals(vertices, faces); 0154 0155 % ensure normals point outward the mesh 0156 if meshVolume(vertices, faces) < 0 0157 normals = -normals; 0158 end 0159 0160 % inner product of normals 0161 dp = sum(normals(ef1, :) .* normals(ef2, :), 2); 0162 0163 % compute the (unsigned) dihedral angle between the normals of the two 0164 % faces incident to each edge 0165 beta = acos(min(max(dp, -1), 1)); 0166 0167 % relatice orientation of face normals cross product and edge orientation 0168 cp = crossProduct3d(normals(ef1, :), normals(ef2, :)); 0169 si = sign(sum(cp .* edgeVectors, 2)); 0170 0171 % compute signed dihedral angle 0172 beta = beta .* si; 0173 0174 0175 %% Compute tensors 0176 0177 if verbose 0178 disp('compute edge tensors'); 0179 end 0180 0181 % curvature tensor of each edge 0182 T = zeros(3, 3, ne); 0183 for i = 1:3 0184 for j = 1:i 0185 T(i, j, :) = reshape(edgeVectors(:,i) .* edgeVectors(:,j), 1, 1, ne); 0186 T(j, i, :) = T(i, j, :); 0187 end 0188 end 0189 T = bsxfun(@times, T, reshape(d .* beta, [1 1 ne])); 0190 0191 % curvature tensor of each vertex by pooling edge tensors 0192 Tv = zeros(3, 3, nv); 0193 w = zeros(1, 1, nv); 0194 for k = 1:ne 0195 if showProgress 0196 progressbar(k, ne); 0197 end 0198 Tv(:,:,ev1(k)) = Tv(:,:,ev1(k)) + T(:,:,k); 0199 Tv(:,:,ev2(k)) = Tv(:,:,ev2(k)) + T(:,:,k); 0200 w(:,:,ev1(k)) = w(:,:,ev1(k)) + 1; 0201 w(:,:,ev2(k)) = w(:,:,ev2(k)) + 1; 0202 end 0203 w(w < eps) = 1; 0204 Tv = Tv ./ repmat(w, [3 3 1]); 0205 0206 if verbose 0207 disp('average vertex tensors'); 0208 end 0209 0210 % apply smoothing on the tensor field 0211 for i = 1:3 0212 for j = 1:3 0213 a = Tv(i, j, :); 0214 a = smoothMeshFunction(vertices, faces, a(:), nIters); 0215 Tv(i, j, :) = reshape(a, [1 1 nv]); 0216 end 0217 end 0218 0219 0220 %% Retrieve curvatures and eigen vectors from tensors 0221 0222 if verbose 0223 disp('retrieve curvatures'); 0224 end 0225 0226 % allocate memory 0227 U = zeros(3, 3, nv); 0228 D = zeros(3, nv); 0229 0230 % iterate over vertices 0231 for k = 1:nv 0232 % display progress 0233 if showProgress 0234 progressbar(k,nv); 0235 end 0236 0237 % extract eigenvectors and eigenvalues for current vertex 0238 [u, d] = eig(Tv(:,:,k)); 0239 d = real(diag(d)); 0240 0241 % sort acording to [normal, min curv, max curv] 0242 [~, I] = sort(abs(d)); 0243 D(:, k) = d(I); 0244 U(:, :, k) = real(u(:,I)); 0245 end 0246 0247 % retrieve main curvatures and associated directions 0248 C1 = D(2,:)'; 0249 C2 = D(3,:)'; 0250 U1 = squeeze(U(:,3,:))'; 0251 U2 = squeeze(U(:,2,:))'; 0252 0253 % enforce C1 < C2 0254 inds = find(C1 > C2); 0255 C1tmp = C1; 0256 U1tmp = U1; 0257 C1(inds) = C2(inds); 0258 C2(inds) = C1tmp(inds); 0259 U1(inds,:) = U2(inds,:); 0260 U2(inds,:) = U1tmp(inds,:); 0261 0262 % compute optional output arguments 0263 if nargout > 4 0264 % average and gaussian curvatures 0265 H = (C1 + C2) / 2; 0266 K = C1 .* C2; 0267 0268 if nargout > 6 0269 % normal vector for each vertex 0270 N = squeeze(U(:,1,:))'; 0271 end 0272 end