面积加权插值python

 

三角网格模型

三角网格模型M通常可用一对线性表表示，M={V,F},V={v;1≤1≤i≤nv}是其顶点表，F={fk:1≤k≤nF}是三角面片表

图1表示顶点

vi的一邻域。图中除顶点 vi外的其他顶点组成的集合记为 Vi。若顶点 vj∈Vi，则 vj是 vi的相邻点。 Vi中的顶点个数记为 |Vi|。包含 vi的三角片集合记为 Fi。若三角片 fk∈Fi，则 fi与 vi是相关联的。 Fi中的三角片个数记为 |Fi|。

三角网格模型顶点法矢计算

图1中三角片fk的法矢Nfk可用下式计算

Nfk=ei,j+1×ei,j||ei,j+1×ei,j||=(vi−vj+1)×(vi−vj)||(vi−vj+1)×(vi−vj)||

式中

ei,j+1表示由顶点 vj指向顶点 vi的边矢量，如图1所示。

计算顶点

vi的法矢 Nvi时，现有文献中常用的方法是将其一邻域内三角片的法矢进行面积加权平均，即按如下公式进行计算：

Nvi=∑fk∈FiAfkNfk||∑fk∈FiAfkNfk||

式中，

Afk表示三角片 fk的面积。

以上算法在我

之前的博文中实现过

如图2所示，当与顶点

vi临接的两个三角片 fk1与 fk2具有相同的面积与法矢时，若用上式计算， fk1与 fk2具有相同的贡献。然而，由图2可知， fk1与 fk2的形状相差很大， p1、 p2分别是两个三角片中离 vi最远的点，然而 p2与 vi的距离要远远大于 p1与 vi的距离。可见，上式未能反映出三角片形状对 Nvi的影响，需要对其进行改进

我们对上式做如下的修正，来计算顶点

vi处的法矢 Nvi：

Nvi=∑fk∈FiγkAfkNfk||∑fk∈FiγkAfkNfk||

式中，

γk为三角片 fk在顶点 vi处的内角，如图1及图2中所示。

因此上式是面积与顶角加权的顶点法矢的计算公式，反映了三角面片面积与形状的综合影响。

MATLAB实现

实现代码如下

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

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

nface = size(face,2);
nvert = size(vertex,2);

% unit normals to the faces
normalf = crossp( vertex(:,face(2,:))-vertex(:,face(1,:)), ...
                  vertex(:,face(3,:))-vertex(:,face(1,:)) );

% unit normal to the vertex
normal = zeros(3,nvert);
for i=1:nface
    f = face(:,i);
    theta=zeros(1,3);
    theta(1)=angle(vertex(:,f(2))-vertex(:,f(1)),vertex(:,f(3))-vertex(:,f(1)));
    theta(2)=angle(vertex(:,f(1))-vertex(:,f(2)),vertex(:,f(3))-vertex(:,f(2)));
    theta(3)=angle(vertex(:,f(1))-vertex(:,f(3)),vertex(:,f(2))-vertex(:,f(3)));
    for j=1:3
        normal(:,f(j)) = normal(:,f(j)) + theta(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
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 theta =angle(a,b)
c=dot(a,b);

d=dot(a,a);  %求a的长度  
e=sqrt(d);  

f=dot(b,b);  %求b的长度  
g=sqrt(f);  

h=c./(e.*g);  

theta=acos(h);
end