网格离散曲率算法（利用Normal cycle 理论计算）

文章理论来自于

Restricted Delaunay Triangulations and Normal Cycle——by David Cohen-Steiner, Jean-Marie Morvan

理论基础

大致是根据利用Normal cycle理论计算曲面的第二基本形式，在给定顶点的周围选取一个小的邻域求平均值。

算法过程

对于网格上的每个顶点$p_i$，取点$p_i$周围的一个测地邻域B,其取法为满足所有到点$p_i$的测地距离小于等于给定距离$r_i$的点，为了简化，可以直接取以点$p_i$为球心，半径为$r_i$的球内所有的点，下面只考虑球内的点和边，定义如下矩阵

$$E_{p_i}(B)=frac{1}{|B|}sum_{ein B}beta(e)||ecap B||bar{e}^{T}cdotbar{e}$$
其中 $|B|$为B的面积， $e$为B中网格的边，$bar{e}$为 $e$方向上的单位向量，$||ecap B||$是 $ecap B$的长度，总是在0和 $|e|$之间， $beta(e)$是以 $e$为公共边的两个法向三角形法向的夹角。如下图所示

考察上面矩阵的特征值和特征方向，$E_{p_i}(B)$的最小特征值和特征方向可作为顶点 $p_i$处网格曲面的法向量的估计式。最小特征值和最大特征值可作为主曲率 $k_1$, $k_2$的估计式，于是平均曲率和高斯曲率就可以直接由主曲率依照微分几何的理论来求得
$$k_G=k_1cdot k_2, k_H=frac{k_1+k_2}{2}$$

代码实现

代码主要还是来自Toolbox graph工具箱，我大致读了一下这段代码，应该只是选择了与目标点直接相连接的点作为邻域，具体代码如下

function [Umin,Umax,Cmin,Cmax,Cmean,Cgauss,Normal] = compute_curvature(vertex,face)

% compute_curvature - compute principal curvature directions and values
%
%   [Umin,Umax,Cmin,Cmax,Cmean,Cgauss,Normal] = compute_curvature(vertex,face,options);
%
%   Umin is the direction of minimum curvature
%   Umax is the direction of maximum curvature
%   Cmin is the minimum curvature
%   Cmax is the maximum curvature
%   Cmean=(Cmin+Cmax)/2
%   Cgauss=Cmin*Cmax
%   Normal is the normal to the surface
%
%   options.curvature_smoothing controls the size of the ring used for
%       averaging the curvature tensor.
%
%   The algorithm is detailed in 
%       David Cohen-Steiner and Jean-Marie Morvan. 
%       Restricted Delaunay triangulations and normal cycle. 
%       In Proc. 19th Annual ACM Symposium on Computational Geometry, 
%       pages 237-246, 2003. 
%   and also in
%       Pierre Alliez, David Cohen-Steiner, Olivier Devillers, Bruno Le巚y, and Mathieu Desbrun. 
%       Anisotropic Polygonal Remeshing. 
%       ACM Transactions on Graphics, 2003. 
%       Note: SIGGRAPH '2003 Conference Proceedings
%
%   Copyright (c) 2007 Gabriel Peyre

orient = 1;


[vertex,face] = check_face_vertex(vertex,face);

n = size(vertex,2);
m = size(face,2);

% associate each edge to a pair of faces 
i = [face(1,:) face(2,:) face(3,:)]; 
j = [face(2,:) face(3,:) face(1,:)];
L = sub2ind([n, n], i, j); 
[a,b]=hist(L,unique(L)); 
Q = ismember(L,b(a>1)); 
[i1, j1] = ind2sub([n, n], L(Q)); 
W = [i' j']; 
R = ismember(W, [i1' j1'], 'rows'); 
J = find(R); 
[Q, K] = sortrows(W(R,:)); 
J = J(K); 
[~, ia, ~] = unique(Q, 'rows'); 
R = true(size(J)); 
R(ia) = false; 
J = J(R);
s = [1:m 1:m 1:m]; 
s(J) = 0; 
A = sparse(i,j,s,n,n); % if (i,j) are an edge at some face then A(i,j) is 
                    % this face index.
[i,j,~] = find(A); % direct link
X = [i j]; 
Y = [j i]; 
[~, I1, ~]=intersect(X, Y, 'rows'); 
I = X(I1, :); 
K1 = I(I(:,1) < I(:,2),:); % only directed edges 
K2 = [K1(:,2) K1(:,1)]; 
[~, ~, s1]= find(A(sub2ind(size(A),K1(:,1), K1(:,2)))); 
[~, ~, s2]= find(A(sub2ind(size(A),K2(:,1), K2(:,2))));
% links edge->faces 
E = [s1 s2]; 
i=K2(:,1); 
j=K2(:,2); 
ne = length(i); % number of directed edges
% normalized edge 
e = vertex(:,j) - vertex(:,i); 
d = sqrt(sum(e.^2,1)); 
d(d<eps) = 1; 
e = e ./ repmat(d,3,1); 
% avoid too large numerics 
d = d./mean(d);

% normals to faces
[~,normal] = compute_normal(vertex,face);

% inner product of normals
dp = sum( normal(:,E(:,1)) .* normal(:,E(:,2)), 1 );
% angle un-signed
beta = acos(clamp(dp,-1,1));
% sign
cp = crossp( normal(:,E(:,1))', normal(:,E(:,2))' )';
si = orient * sign( sum( cp.*e,1 ) );
% angle signed
beta = beta .* si;
% tensors
T = zeros(3,3,ne);
for x=1:3
    for y=1:x
        T(x,y,:) = reshape( e(x,:).*e(y,:), 1,1,ne );
        T(y,x,:) = T(x,y,:);
    end
end
T = T.*repmat( reshape(d.*beta,1,1,ne), [3,3,1] );

% do pooling on vertices
Tv = zeros(3,3,n);
w = zeros(1,1,n);
for k=1:ne
%    progressbar(k,ne);
    Tv(:,:,i(k)) = Tv(:,:,i(k)) + T(:,:,k);
    Tv(:,:,j(k)) = Tv(:,:,j(k)) + T(:,:,k);
    w(:,:,i(k)) = w(:,:,i(k)) + 1;
    w(:,:,j(k)) = w(:,:,j(k)) + 1;
end
w(w<eps) = 1;
Tv = Tv./repmat(w,[3,3,1]);


% extract eigenvectors and eigenvalues
U = zeros(3,3,n);
D = zeros(3,n);
for k=1:n
    [u,d] = eig(Tv(:,:,k));
    d = real(diag(d));
    % sort acording to [norma,min curv, max curv]
    [~,I] = sort(abs(d));    
    D(:,k) = d(I);
    U(:,:,k) = real(u(:,I));
end

Umin = squeeze(U(:,3,:));
Umax = squeeze(U(:,2,:));
Cmin = D(2,:)';
Cmax = D(3,:)';
Normal = squeeze(U(:,1,:));
Cmean = (Cmin+Cmax)/2;
Cgauss = Cmin.*Cmax;

% enforce than min<max
I = find(Cmin>Cmax);
Cmin1 = Cmin; Umin1 = Umin;
Cmin(I) = Cmax(I); Cmax(I) = Cmin1(I);
Umin(:,I) = Umax(:,I); Umax(:,I) = Umin1(:,I);

% try to re-orient the normals
normal = compute_normal(vertex,face);
s = sign( sum(Normal.*normal,1) ); 
Normal = Normal .* repmat(s, 3,1);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function z = crossp(x,y)
% x and y are (m,3) dimensional
z = x;
z(:,1) = x(:,2).*y(:,3) - x(:,3).*y(:,2);
z(:,2) = x(:,3).*y(:,1) - x(:,1).*y(:,3);
z(:,3) = x(:,1).*y(:,2) - x(:,2).*y(:,1);
end

function y = clamp(x,a,b)

% clamp - clamp a value
%
%   y = clamp(x,a,b);
%
% Default is [a,b]=[0,1].
%
%   Copyright (c) 2004 Gabriel Peyr

if nargin<2
    a = 0;
end
if nargin<3
    b = 1;
end

y = max(x,a);
y = min(y,b);
end

其中check_face_vertex.m文件为

function [vertex,face] = check_face_vertex(vertex,face)

% check_face_vertex - check that vertices and faces have the correct size
%
%   [vertex,face] = check_face_vertex(vertex,face);
%
%   Copyright (c) 2007 Gabriel Peyre

vertex = check_size(vertex,2,4);
face = check_size(face,3,4);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function a = check_size(a,vmin,vmax)
if isempty(a)
    return;
end
if size(a,1)>size(a,2)
    a = a';
end
if size(a,1)<3 && size(a,2)==3
    a = a';
end
if size(a,1)<=3 && size(a,2)>=3 && sum(abs(a(:,3)))==0
    % for flat triangles
    a = a';
end
if size(a,1)<vmin ||  size(a,1)>vmax
    error('face or vertex is not of correct size');
end
end

还有compute_normal.m文件计算平面法向

function [normal,normalf] = compute_normal(vertex,face)

% compute_normal - compute the normal of a triangulation
%
%   [normal,normalf] = compute_normal(vertex,face);
%
%   normal(i,:) is the normal at vertex i.
%   normalf(j,:) is the normal at face j.
%
%   Copyright (c) 2004 Gabriel Peyr

[vertex,face] = check_face_vertex(vertex,face);

nface = size(face,2);
nvert = size(vertex,2);
normal = zeros(3, nvert);

% unit normals to the faces
normalf = crossp( vertex(:,face(2,:))-vertex(:,face(1,:)), ...
                  vertex(:,face(3,:))-vertex(:,face(1,:)) );
d = sqrt( sum(normalf.^2,1) ); d(d<eps)=1;
normalf = normalf ./ repmat( d, 3,1 );

% unit normal to the vertex
normal = zeros(3,nvert);
for i=1:nface
    f = face(:,i);
    for j=1:3
        normal(:,f(j)) = normal(:,f(j)) + normalf(:,i);
    end
end
% normalize
d = sqrt( sum(normal.^2,1) ); d(d<eps)=1;
normal = normal ./ repmat( d, 3,1 );

% enforce that the normal are outward
v = vertex - repmat(mean(vertex,1), 3,1);
s = sum( v.*normal, 2 );
if sum(s>0)<sum(s<0)
    % flip
    normal = -normal;
    normalf = -normalf;
end



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function z = crossp(x,y)
% x and y are (m,3) dimensional
z = x;
z(1,:) = x(2,:).*y(3,:) - x(3,:).*y(2,:);
z(2,:) = x(3,:).*y(1,:) - x(1,:).*y(3,:);
z(3,:) = x(1,:).*y(2,:) - x(2,:).*y(1,:);

演示

在以小兔子为例，首先展示高斯曲率和平均曲率，然后就是一段执行代码了

[vertices,faces]=obj__read('bunny_200.obj');
[Umin,Umax,Cmin,Cmax,Cmean,Cgauss,Normal] = compute_curvature(vertices,faces);
trep=triangulation(faces',vertices');

绘制高斯曲率图

colormap jet;trisurf(trep,Cgauss);shading interp;caxis([min(Cgauss) max(Cgauss)]);axis square;colorbar('vert');brighten(-0.1);axis off;

绘制平均曲率图

colormap jet;trisurf(trep,Cmean);shading interp;caxis([min(Cmean) max(Cmean)]);axis square;colorbar('vert');brighten(-0.1);axis off;

并且根据测试，所有的值是不依赖于坐标系的（内蕴的）
再对两个主法向进行绘制

figure, title('Principal A');
p1=vertices-0.02*Umin; p2=vertices+0.02*Umin;
p1=p1';p2=p2';
plot3([p1(:,1) p2(:,1)]',[p1(:,2) p2(:,2)]',[p1(:,3) p2(:,3)]','k-');hold on;trimesh(trep);axis off;
axis equal; view(3);

figure, title('Principal B');
p1=vertices-0.02*Umax; p2=vertices+0.02*Umax;
p1=p1';p2=p2';
plot3([p1(:,1) p2(:,1)]',[p1(:,2) p2(:,2)]',[p1(:,3) p2(:,3)]','k-');hold on;trimesh(trep);axis off;
axis equal; view(3);

再做一个小例子，人脸

colormap jet;trisurf(trep,Cgauss);shading interp;caxis([min(Cgauss) max(Cgauss)-0.1]);axis square;colorbar('vert');brighten(-0.1);axis off;

colormap jet;trisurf(trep,Cmean);shading interp;caxis([min(Cmean) max(Cmean)]);axis square;colorbar('vert');brighten(-0.1);axis off;

分别对两幅图的上限减少了0.1（用来使得变化快的部分看起来更明显）

gauss曲率图	平均曲率图

最后

以上就是缓慢芒果为你收集整理的网格离散曲率算法（利用Normal cycle 理论计算）的全部内容，希望文章能够帮你解决网格离散曲率算法（利用Normal cycle 理论计算）所遇到的程序开发问题。

如果觉得靠谱客网站的内容还不错，欢迎将靠谱客网站推荐给程序员好友。