对于网格曲面,顶点法矢计算的常用方法是将其一邻域内三角片的法矢进行面积加权平均。
这里参考下面文献

三角网格模型顶点法矢与离散曲率计算 ——神会存,李建华,周来水

提出了一种新的面积加权方法。

三角网格模型

三角网格模型M通常可用一对线性表表示, M={V,F},V={v;11inv} <script type="math/tex" id="MathJax-Element-1">M=\{V,F\},V=\{v;1\leq1\leq i\leq n_v\}</script>是其顶点表, F={fk:1knF} <script type="math/tex" id="MathJax-Element-2">F=\{f_k:1\leq k\leq n_F\}</script>是三角面片表

这里写图片描述

图1表示顶点 vi <script type="math/tex" id="MathJax-Element-3">v_i</script>的一邻域。图中除顶点 vi <script type="math/tex" id="MathJax-Element-4">v_i</script>外的其他顶点组成的集合记为 Vi <script type="math/tex" id="MathJax-Element-5">V^i</script>。若顶点 vjVi <script type="math/tex" id="MathJax-Element-6">v_j\in V^i</script>,则 vj <script type="math/tex" id="MathJax-Element-7">v_j</script>是 vi <script type="math/tex" id="MathJax-Element-8">v_i</script>的相邻点。 Vi <script type="math/tex" id="MathJax-Element-9">V^i</script>中的顶点个数记为 |Vi| <script type="math/tex" id="MathJax-Element-10">|V^i|</script>。包含 vi <script type="math/tex" id="MathJax-Element-11">v_i</script>的三角片集合记为 Fi <script type="math/tex" id="MathJax-Element-12">F^i</script>。若三角片 fkFi <script type="math/tex" id="MathJax-Element-13">f_k\in F^i</script>,则 fi <script type="math/tex" id="MathJax-Element-14">f_i</script>与 vi <script type="math/tex" id="MathJax-Element-15">v_i</script>是相关联的。 Fi <script type="math/tex" id="MathJax-Element-16">F^i</script>中的三角片个数记为 |Fi| <script type="math/tex" id="MathJax-Element-17">|F^i|</script>。

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

图1中三角片 fk <script type="math/tex" id="MathJax-Element-18">f_k</script>的法矢 Nfk <script type="math/tex" id="MathJax-Element-19">N_{f_k}</script>可用下式计算

Nfk=ei,j+1×ei,j||ei,j+1×ei,j||=(vivj+1)×(vivj)||(vivj+1)×(vivj)||
<script type="math/tex; mode=display" id="MathJax-Element-20">N_{f_k}=\frac{e_{i,j+1}\times e_{i,j}}{||e_{i,j+1}\times e_{i,j}||}=\frac{(v_i-v_{j+1})\times(v_i-v_j)}{||(v_i-v_{j+1})\times(v_i-v_j)||}</script>
式中 ei,j+1 <script type="math/tex" id="MathJax-Element-21">e_{i,j+1}</script>表示由顶点 vj <script type="math/tex" id="MathJax-Element-22">v_j</script>指向顶点 vi <script type="math/tex" id="MathJax-Element-23">v_i</script>的边矢量,如图1所示。
计算顶点 vi <script type="math/tex" id="MathJax-Element-24">v_i</script>的法矢 Nvi <script type="math/tex" id="MathJax-Element-25">N_{v_i}</script>时,现有文献中常用的方法是将其一邻域内三角片的法矢进行面积加权平均,即按如下公式进行计算:
Nvi=fkFiAfkNfk||fkFiAfkNfk||
<script type="math/tex; mode=display" id="MathJax-Element-26">N_{v_i}=\frac{\sum_{f_k \in F^i}A_{f_k}N_{f_k}}{||\sum_{f_k \in F^i}A_{f_k}N_{f_k}||}</script>
式中, Afk <script type="math/tex" id="MathJax-Element-27">A_{f_k}</script>表示三角片 fk <script type="math/tex" id="MathJax-Element-28">f_k</script>的面积。
以上算法在我 之前的博文中实现过

这里写图片描述

如图2所示,当与顶点 vi <script type="math/tex" id="MathJax-Element-29">v_i</script>临接的两个三角片 fk1 <script type="math/tex" id="MathJax-Element-30">f_{k_1}</script>与 fk2 <script type="math/tex" id="MathJax-Element-31">f_{k_2}</script>具有相同的面积与法矢时,若用上式计算, fk1 <script type="math/tex" id="MathJax-Element-32">f_{k_1}</script>与 fk2 <script type="math/tex" id="MathJax-Element-33">f_{k_2}</script>具有相同的贡献。然而,由图2可知, fk1 <script type="math/tex" id="MathJax-Element-34">f_{k_1}</script>与 fk2 <script type="math/tex" id="MathJax-Element-35">f_{k_2}</script>的形状相差很大, p1 <script type="math/tex" id="MathJax-Element-36">p_1</script>、 p2 <script type="math/tex" id="MathJax-Element-37">p_2</script>分别是两个三角片中离 vi <script type="math/tex" id="MathJax-Element-38">v_i</script>最远的点,然而 p2 <script type="math/tex" id="MathJax-Element-39">p_2</script>与 vi <script type="math/tex" id="MathJax-Element-40">v_i</script>的距离要远远大于 p1 <script type="math/tex" id="MathJax-Element-41">p_1</script>与 vi <script type="math/tex" id="MathJax-Element-42">v_i</script>的距离。可见,上式未能反映出三角片形状对 Nvi <script type="math/tex" id="MathJax-Element-43">N_{v_i}</script>的影响,需要对其进行改进
我们对上式做如下的修正,来计算顶点 vi <script type="math/tex" id="MathJax-Element-44">v_i</script>处的法矢 Nvi <script type="math/tex" id="MathJax-Element-45">N_{v_i}</script>:
Nvi=fkFiγkAfkNfk||fkFiγkAfkNfk||
<script type="math/tex; mode=display" id="MathJax-Element-46">N_{v_i}=\frac{\sum_{f_k \in F^i}\gamma_k A_{f_k}N_{f_k}}{||\sum_{f_k \in F^i}\gamma_kA_{f_k}N_{f_k}||}</script>
式中, γk <script type="math/tex" id="MathJax-Element-47">\gamma_k</script>为三角片 fk <script type="math/tex" id="MathJax-Element-48">f_k</script>在顶点 vi <script type="math/tex" id="MathJax-Element-49">v_i</script>处的内角,如图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

完成~

Logo

DAMO开发者矩阵,由阿里巴巴达摩院和中国互联网协会联合发起,致力于探讨最前沿的技术趋势与应用成果,搭建高质量的交流与分享平台,推动技术创新与产业应用链接,围绕“人工智能与新型计算”构建开放共享的开发者生态。

更多推荐