一种网格去噪算法（基于平均面法向的均值滤波）
程序员文章站
2024-03-25 22:27:58
...
算法原文来自
Mesh smoothing via mean and median filtering applied to face normals——H. Yagou, Y. Ohtake, and A. Belyaev
MathWorks论坛中有Two functions for smoothing/denoising of triangular meshes给出了这种算法的Matlab代码，下面介绍一下算法的步骤
符号说明及算法步骤
符号说明
算法步骤
- 对于每个网格三角形
ti ,计算区域加权的面法向m(ti) m(ti)=1∑tj∈t∗iA(tj)∑tj∈t∗iA(tj)n(tj) - 规范化平均法向
m(ti) m(ti)←m(ti)||m(ti)|| - 用如下方法更新网格中每个点的坐标
vi←vi+1∑tj∈t∗iA(tj)∑tj∈t∗iA(tj)πi(tj)
其中πi(tj)=<eij,m(tj)>m(tj) 并且eij=cj−vi 是从顶点vi 到三角形tj 中心cj 的向量。注意到由内积的定义，向量πi(tj) 是向量eij 在法向tj 上的投影。
Matlab代码
meannorm_trismooth.m函数
function xyzn=meannorm_trismooth(xyz,t)
% Mean face normal filter for smoothing/denoising triangular meshes
%
% Reference: 1) Yagou, Belayev, Ohtake (2002) Mesh smoothing via Mean and
% Median Filtering applied to face Normals, PGMP - Theory and
% Applications
% 2) Zhang and Hamza, (2006) Vertex-based anisotropic smoothing of
% 3D mesh data, IEEE CCECE
% Acknowledgement:
% Q. Fang: iso2mesh (http://iso2mesh.sf.net)
%
% Input: xyz <nx3> vertex coordinates
% t <mx3> triangulation index array
% Output: xyzn <nx3> updated vertex coordinates
% Version: 1
% JOK 300709
% I/O check:
if nargin~=2
error('Wrong # of input')
end
if nargout ~= 1
error('Output is a single array, wrong designation!')
end
nt= size(t);
if nt(2)~=3
error('Triangle element matrix should be mx3!')
end
mp = size(xyz);
if mp(2)~=3
error('Vertices should be nx3!')
end
% Containing the elements adjacent to each vertex within t
[conn,connnum,count]=neighborelem(t,max(max(t)));% Triangles adjacent to each vertex;
trineigh = trineighborelem(t);% Containing each element adjacent to each elements
tarea = triangle_area(xyz,t);
% Initialize etc.
xyzn = zeros(size(xyz));
nvec = trinormal(t,xyz);
nvecn = zeros(size(nvec));
tol = 1e-2;err = Inf;
itermax = 30;iter = 0;
while err>tol && iter<itermax
iter = iter+1;
for k = 1:length(t)
% Find triangles neighborhood
indt01=trineigh{k}; % Element indices
% Step 1: Area weighted face normal
mti=1/sum(tarea(indt01))*sum(repmat(tarea(indt01),1,3).*nvec(indt01,:),1);
% Step 2: Normalize mti & update
lmti = sqrt(sum(mti.*mti,2));
mti = mti./repmat(lmti,1,3);
nvecn(k,:) = mti;
% Step 3: update each vertex
if length(indt01)>2
for i = 1:3
% Evaluate for each vertex in the central triangle: t(k,:)
indvaux01 = t(k,i); % Current vertex
indvaux02=conn{indvaux01};% Elements containing current vertex
cj = tricentroid(xyz,t(indvaux02,:));
nc=size(cj);
eij = cj-repmat(xyz(t(k,i),:),nc(1),1);
pii = repmat(dot(eij,nvecn(indvaux02,:),2),1,3).*nvecn(indvaux02,:);
xyzn(t(k,i),:)=xyz(t(k,i),:)+1/sum(tarea(indvaux02))*sum(repmat(tarea(indvaux02),1,3).*pii,1);
end% for
else
xyzn(t(k,:),:) = xyzn(t(k,:),:);
end% if
end% for k
tarean = triangle_area(xyzn,t);
% Stopping criteria: face-normal error metric (L2-norm regardless of
% median or mean filter used above)
err = 1/sum(tarean)*sum(tarean.*sqrt(sum((nvec-nvecn).*(nvec-nvecn),2)),1);
% Update coordinates, triangle normal vectors, mesh areas
xyz = xyzn;
nvec = nvecn;
tarea=tarean;
end%while iter/err
end % meannormfilter_smooth_v01
% %%%%%%%%%%%%%%%%%%%%%% SUBFUNCTIONS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
function [conn,connnum,count]=neighborelem(elem,nn)
% [conn,connnum,count]=neighborelem(elem,nn)
% neighborelem: create node neighbor list from a mesh
% parameters:
% elem: element table of a mesh
% nn : total node number of the mesh
% conn: output, a cell structure of length nn, conn{n}
% contains a list of all neighboring elem ID for node n
% connnum: vector of length nn, denotes the neighbor number of each node
% count: total neighbor numbers
conn=cell(nn,1);
dim=size(elem);
for i=1:dim(1)
for j=1:dim(2)
conn{elem(i,j)}=[conn{elem(i,j)},i];
end
end
count=0;
connnum=zeros(1,nn);
for i=1:nn
conn{i}=sort(conn{i});
connnum(i)=length(conn{i});
count=count+connnum(i);
end
end % neighborelem
%%
function [area]=triangle_area(xyz,t)
% This function gives the area of a triangle
aux1 = xyz(t(:,3),:)-xyz(t(:,1),:);
aux2 = xyz(t(:,2),:)-xyz(t(:,1),:);
aux3 = xyz(t(:,3),:)-xyz(t(:,2),:);
L = sort([sqrt(sum(aux1.*aux1,2)),sqrt(sum(aux2.*aux2,2)),sqrt(sum(aux3.*aux3,2))],2);
% Heron's formula (stable)
area= sqrt( (L(:,3)+(L(:,2)+L(:,1))).*(L(:,3)+(L(:,2)-L(:,1))).*(L(:,1)+(L(:,3)-L(:,2))).*(L(:,1)-(L(:,3)-L(:,2))) )./4;
end % triangle_area
%%
function out = tricentroid(v,tri)
% Function to output the centroid of triangluar elements.
% Note that the output will be of length(tri)x3
% Input: <v> nx2 or 3: vertices referenced in tri
% <tri> mx3: triangle indices
% Version: 1
% JOK 300509
% I/O check
[nv,mv]=size(v);
[nt,mt]=size(tri);
if mv==2
v(:,3) = zeros(nv,1);
elseif mt~=3
tri=tri';
end
out(:,1) = 1/3*(v(tri(:,1),1)+v(tri(:,2),1)+v(tri(:,3),1));
out(:,2) = 1/3*(v(tri(:,1),2)+v(tri(:,2),2)+v(tri(:,3),2));
out(:,3) = 1/3*(v(tri(:,1),3)+v(tri(:,2),3)+v(tri(:,3),3));
end% tricentroid
%%
function trineigh=trineighborelem(elem)
% [conn,connnum,count]=neighborelem(elem,nn)
%
% neighborelem: create node neighbor list from a mesh
%
% author: fangq (fangq<at> nmr.mgh.harvard.edu)
% date: 2007/11/21
%
% parameters:
% elem: element table of a mesh
% nn : total node number of the mesh
% conn: output, a cell structure of length nn, conn{n}
% contains a list of all neighboring elem ID for node n
% connnum: vector of length nn, denotes the neighbor number of each node
% count: total neighbor numbers
%
% -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)
%
nn = max(max(elem));
conn=cell(nn,1);
dim=size(elem);
for i=1:dim(1)
for j=1:dim(2)
conn{elem(i,j)}=[conn{elem(i,j)},i];
end
end
count=0;
connnum=zeros(1,nn);
for i=1:nn
conn{i}=sort(conn{i});
connnum(i)=length(conn{i});
count=count+connnum(i);
end
% Now identify all the nodes in one triangle and get them all into one
% array, so that all elements surrounding one element are listed
trineigh = cell(dim(1),1);
for i=1:dim(1)
vind = elem(i,:); % Vertices for the i'th element
aux01 = [conn{vind(1)},conn{vind(2)},conn{vind(3)}];
aux02 = unique(aux01);
trineigh{i} = aux02;
end
end%trineighborelem
function nvec=trinormal(tri,p)
% Function to compute normal vector to triangle
% Input: tri <triangular index matrix>
% p nx3 array of vertices
% Output: nvec <nx3> array of normal vector
% JOK180609
% Version: 1
% I/O check
% to be completed
% Check polarity of elements
% Construct vectors
v = [p(tri(:,3),1)-p(tri(:,1),1), p(tri(:,3),2)-p(tri(:,1),2), p(tri(:,3),3)-p(tri(:,1),3)];
w = [p(tri(:,2),1)-p(tri(:,1),1), p(tri(:,2),2)-p(tri(:,1),2), p(tri(:,2),3)-p(tri(:,1),3)];
% Calculate cross product
normvec = [v(:,2).*w(:,3)-v(:,3).*w(:,2), ...
-(v(:,1).*w(:,3)-v(:,3).*w(:,1)), ...
v(:,1).*w(:,2)-v(:,2).*w(:,1)];
% Normalize
lnvec = sqrt(sum(normvec.*normvec,2));
nvec = normvec./repmat(lnvec,1,3);
end演示代码及实例
下面展示了龙以及误差分布效果
function test
%TEST Summary of this function goes here
% Detailed explanation goes here
[vertex,face]=obj__read('Dragon2002.obj');
normals = compute_normal(vertex,face);
rho = randn(1,size(vertex,2))*0.03;
vertex1 = vertex + repmat(rho,3,1).*normals;
vertex2 = vertex1';
for i=1:3
subplot(1,3,i);
trep=triangulation(face',vertex2);
vertex2=meannorm_trismooth(vertex2,face');
% %trimesh(trep);axis square;axis off;
d=sqrt(sum((vertex2-vertex').^2,2));
colormap jet;trisurf(trep,d);shading interp;caxis([0 0.05]);axis square;colorbar('vert');brighten(-0.1);axis off;
end
end下面展示误差分布
| 加噪声的龙 | 迭代一次 | 迭代两次 |
如果是绘制网格的话，上面的语句改为
for i=1:3
subplot(1,3,i);
trep=triangulation(face',vertex2);
vertex2=meannorm_trismooth(vertex2,face');
trimesh(trep);axis square;axis off;下面展示误差分布
| 加噪声的龙 | 迭代一次 | 迭代两次 |
个人感觉效果并不如之前的算法（可能多试几个例子就好了？）