Libigl 包含各种各样的几何处理工具和函数,用于处理网格以及与其相关的线性代数:在本介绍性教程中讨论的内容太多了。 我们在本章中提取了一些有趣的函数来强调。
网格统计
Libigl 包含各种网格统计信息,包括面角度、面面积和奇异顶点的检测,奇异顶点是三角剖分中具有多于或少于 6 个邻居或四边形中具有多于或少于 4 个邻居的顶点。
示例统计计算这些数量并进行基本的统计分析,以估计网格的等轴度和规则性:
Irregular vertices:
136/2400 (5.67%)
Areas (Min/Max)/Avg_Area Sigma:
0.01/5.33 (0.87)
Angles in degrees (Min/Max) Sigma:
17.21/171.79 (15.36)
第一行包含不规则顶点的数量和百分比,这对于四边形网格用于定义细分曲面时尤为重要:每个奇异点将导致曲面上仅存在一个点
c
∧
1
c^\wedge1
c∧1.
第二行报告最小元素、最大元素的面积和标准差。 这些数字通过平均面积进行归一化,因此在上面的示例中,5.33 最大面积意味着最大的面片比平均面片大 5 倍。 理想的各向同性网格的最小和最大面积都接近 1。
第三行测量面角,在完美规则的三角测量中,面角应接近 60 度(四边形为 90 度)。 对于 FEM(finite element method) 目的,角度越接近 60 度,优化就越稳定。 在这种情况下,很明显网格质量很差,如果用于求解偏微分方程,可能会导致错误。
广义绕组数
封闭水密表面网格内部的四面体化问题是一个困难但适定的问题(请参阅我们的 Tetgen 包装器)。 但是像 TetGen 这样的黑盒四面体网格划分器将拒绝输入具有自相交、开放边界、来自多个连接组件的非流形边缘的三角形网格。 问题有两个方面:自相交提出了矛盾的面约束,并且自相交/开放边界/非流形边缘使得在没有进一步假设的情况下从外部确定内部的问题变得不适定。
第一个问题很容易通过“解决”所有自相交来解决。 也就是说,对相交的三角形进行网格划分,以便相交恰好发生在边和顶点处。 这是使用 igl::selfintersect 完成的。
TetGen 通常可以对这个“已解析”网格的凸包进行四面体化,然后问题就变成确定这些四面体中哪些在输入网格内部,哪些在外部。 也就是说,哪些应该保留,哪些应该删除。
“广义缠绕数”是一种用于确定麻烦网格的内部和外部的可靠方法
在某个点
p
∈
R
3
p\in R^3
p∈R3关于(V,F)的广义绕数定义为标量函数:
ω
(
p
)
=
∑
f
i
∈
F
1
4
π
Ω
f
i
(
p
)
\omega(p)=\sum_{f_i\in F}\frac{1}{4\pi}\Omega_{f_i}(p)
ω(p)=∑fi∈F4π1Ωfi(p)
其中
Ω
f
i
\Omega_{f_i}
Ωfi是
f
i
f_i
fi(F 中的第 i 个面)在该点p所对的立体角。该立体角贡献是一个简单的封闭式表达式,涉及 atan2 和一些点积。
如果(V,F)是一个封闭的水密表面,那么如果点p在(V,F)内部,有
ω
(
p
)
=
1
\omega(p)=1
ω(p)=1,外部则有
ω
(
p
)
=
0
\omega(p)=0
ω(p)=0。如果 (V,F) 闭合但自身重叠,那么如果 (V,F) 闭合但自身重叠
ω
(
p
)
\omega(p)
ω(p)是一个整数,这个整数代表着(V,F)包围顶点p的次数。如果(V,F)不是一个闭合的表面或者根本不是一个流形的表面(但是它的法向起码是一致的),那么如果点p接近于(V,F)的内部,那么
ω
(
p
)
\omega(p)
ω(p)趋近于1,如果点p接近于(V,F)的外部,那么
ω
(
p
)
\omega(p)
ω(p)趋近于0。
示例 702 计算具有孔和自交点的猫内部四面体网格的广义缠绕数函数(金色)。 银色网格是提取的内部四面体的表面,切片显示凸包中所有四面体上的缠绕数函数:蓝色(〜0),绿色(〜1),黄色(〜2)。
网格抽取(简化)
网格简化或抽取的研究几乎与网格本身一样古老。 给定一个具有太多三角形的高分辨率网格,找到一个“非常近似”的具有更少三角形的低分辨率网格。 到目前为止,有多种不同的范式可以解决这个问题,而且最先进的方法也相当先进。
一种网格抽取方法通过从网格中连续移除元素来进行操作。 特别是,Hoppe 提倡连续删除或折叠边。 该技术的通用形式是通过折叠一条边构造一个由初始高分辨率网格
M
0
M_0
M0 到最低分辨率网格
M
n
M_n
Mn 组成的 n 个网格序列:
M
0
→
e
d
g
e
c
o
l
l
a
p
s
e
M
1
→
e
d
g
e
c
o
l
l
a
p
s
e
.
.
.
→
e
d
g
e
c
o
o
l
a
p
s
e
M
n
−
1
→
e
d
g
e
c
o
l
l
a
p
s
e
M
n
M_0\to_{edge\ collapse} M_1\to_{edge\ collapse}...\to_{edge\ coolapse}M_{n-1}\to_{edge\ collapse}M_n
M0→edge collapseM1→edge collapse...→edge coolapseMn−1→edge collapseMn
Hoppe 的原始方法和随后的后续工作提出了各种方法来选择此序列中下一个要折叠的边。 使用基于成本的范例,可以根据边缘的“成本”维护边缘的优先级队列(如果删除该边缘,我的近似值会“更差”多少?)。 最便宜的边被折叠,相邻边的成本被更新。
为了维持拓扑(例如,如果网格组合为球体或环面等),应该为边缘分配无限成本,这些边缘的折叠会改变网格拓扑。 事实上,当且仅当折叠边端点的相互邻居的数量不正好是两个时,才会发生这种情况!
如果存在第三个共享顶点,则将删除另一个面,而且删除 了2 条边。 这可能会导致出现不需要的孔或非流形的翻转。
还有一个一次性条件,即四面体的任何边都不应塌陷。
由于 libigl(故意)不以动态网格数据结构(例如半边数据结构)为中心实现其实现,因此对拓扑更改的支持受到限制。 尽管如此,libigl 支持孤立的边缘塌陷、边缘塌陷序列(每个都在 O(log) 时间内)和基于优先级队列的抽取。
最简单的是 igl::decimation。 通过调用
igl::decimate(V,F,1000,U,G);
其中,shortest_edge_and_midpoint 将边的长度指定为成本,将其中点指定为合并的顶点位置,max_m 计算当前的面数(有效折叠将计数减少 2),如果计数低于 m=1000,则返回 true。
我们还可以深入了解抽取循环并直接调用 igl::collapse_edge。 为了高效运行,该例程需要比通常的 (V,F) 网格表示更多的内容。 我们需要 E 一个边索引列表,其中 E.row(i) --> [s,d]; 我们需要 EMAP,它将 F 中每个三角形的“半”边映射到 E 中相应的边,这样 E.row(EMAP(f+i*F.rows)) --> [s,d] 如果边 第 f 个面的第 i 个角的对面是 [s,d](直到方向); 我们需要 EF 和 EI 来跟踪入射在每条边上的面以及边出现在这些面的哪个角上,因此 EF(e,o) = f 和 EI(e,o) = i 意味着该边 E.row(e) --> [s,d] 出现在第 i 个角对面的第 f 个面上(对于 o=0,边缘方向应该匹配,对于 o=1,方向相反)。
当塌陷发生时,F、E等矩阵的大小不会改变。 相反,与“删除”的面和边相对应的行被设置为特殊的常量值 IGL_COLLAPSE_EDGE_NULL。 这样做可以确保我们能够在真正恒定的时间 O(1) 内删除边缘。(以上部分涉及到的拓扑操作为半边数据结构的基本操作)
方便的是 IGL_COLLAPSE_EDGE_NULL==0。 这意味着 F 的大多数 OPENGL 风格渲染将简单地在第一个顶点绘制一堆 0 面积三角形。
以下代码将折叠第一条边并将其合并的顶点放置在原点:
igl::collapse_edge(0,RowVector3d(0,0,0),V,F,E,EMAP,EF,EI);
如果有效(需关注有效性的条件判定),则相应调整 V、F、E、EF、EI。
这个功能很强大,但是级别很低。 要围绕此构建抽取器,您需要跟踪哪些边缘要折叠以及接下来要折叠哪些边缘。 幸运的是,libigl 还公开了基于优先级队列的边缘折叠,并带有函数句柄来调整成本和布局。
优先级队列被实现为(有序)集合 Q 或(成本,边索引)对和迭代器 Qit 列表,以便 Qit[e] 显示 Q 中对应于 eth 边的迭代器。 放置位置存储在位置 C 的 #E 列表中。当调用以下命令时:
igl::collapse_edge(cost_and_placement,V,F,E,EMAP,EF,EI,Q,Qit,C);
尝试根据 Q 的最低成本边缘塌陷。 如果有效,则 V、F 等。 进行相应调整,并且该边会从 Q 中“弹出”。使用 Qit,其相邻边也会从 Q 中弹出,并在根据 cost_and_placement 更新其成本后重新插入,新的位置会被记住在 C 中。如果无效,则该边将被 从 Q 中“弹出”并以无限成本重新插入。
示例 703 演示了使用这种基于优先级队列的方法以及上面讨论的简单的最短边中点成本/布局策略。
符号距离
在“广义缠绕数”部分中,我们研究了一种稳健的方法来确定点是否位于给定三角形汤网格的内部或外部。 Libigl 通过加速有符号和无符号距离查询以及平面三角形网格和 3D 四面体网格的“元素内”查询来补充该算法。 这些例程利用 libigl 的通用轴对齐边界框层次结构 (igl/AABB.h)。 此类是轻量级的,并且根据设计,不存储网格的副本(而是将其作为其成员函数的输入)
顶点位置
对于四面体网格,这对于“元素内”或“点位置”查询很有用:给定一个顶点 q ∈ R 3 q\in R^3 q∈R3以及一个四面体网格 ( V , T ) (V,T) (V,T),判定顶点位于那个四面体当中。在 libigl 中v这是通过 igl::in_element() 针对四面体网格 V,T 和 Q 行中的查询点列表完成的:
// Initialize AABB tree
igl::AABB<MatrixXd,3> tree;
tree.init(V,T);
VectorXi I;
igl::in_element(V,T,Q,tree,I);
得到的向量 I 是 T 中的索引列表,揭示了发现的包含 Q 中对应点的第一个四面体。
对于重叠网格,点 q 可能属于多个四面体。 在这些情况下,可以通过使用 igl::in_element 重载和 SparseMatrix 作为输出来找到所有它们(不仅仅是第一个):
SparseMatrix<int> I;
igl::in_element(V,T,Q,tree,I);
现在 I 的每一行都揭示了每个 tet 是否包含 Q 中的相应行: I(q,e)!=0 表示点 q 在元素 e 中。
最邻近顶点
对于三角形网格,我们使用AABB树来加速点网格最近点查询:给定一个三角网格 ( V , F ) (V,F) (V,F)以及一个待查询的顶点 q ∈ R 3 q\in R^3 q∈R3,找到一个距离q最邻近的顶点 c ∈ ( V , F ) c\in (V,F) c∈(V,F)(其中的顶点c并不一定是三角网格的某个顶点)。这是通过 igl::point_mesh_squared_distance 对三角形网格 V,F 和 P 行中的点列表完成的:
VectorXd sqrD;
VectorXi I;
MatrixXd C;
igl::point_mesh_squared_distance(P,V,F,sqrD,I,C);
输出 sqrD 包含从 P 中的每个点到 C 中给出的最近点的(无符号)平方距离,该点位于由 I 给出的 F 中的元素上(例如,可以使用 igl::barycentric_coordinates 从中恢复重心坐标)。
如果网格 V,F 是静态的,但点集 P 是动态变化的,那么最好重用在 igl::point_mesh_squared_distance 期间构建的 AABB 层次结构:
igl::AABB tree;
tree.init(V,F);
tree.squared_distance(V,F,P,sqrD,I,C);
... // P changes, but (V,F) does not
tree.squared_distance(V,F,P,sqrD,I,C);
符号距离
最后,从最近的点或绕数可以为这个距离施加符号。 在 igl::signed_distance 中,我们提供了两种添加符号方法:所谓的“伪正态测试”37[] 和广义绕数 40[]。
伪法线测试(另请参见 igl::pseudonormal_test)假设输入网格是水密(封闭、非自相交、流形)网格。然后给定一个查询的顶点q以及他的最邻近顶点
c
∈
(
V
,
F
)
c\in (V,F)
c∈(V,F)。仔细地选择 c 处的外向法线 n使得
(
q
−
c
)
⋅
n
(\bold{q-c})\cdot\bold n
(q−c)⋅n的符号表示点q是否在三角网格的内部,如果这个符号是+1,顶点在网格外部,若是-1,顶点在网格的内部。如果顶点c的位置确定了,这是一个O(1)的快速测试,但是这个测试会在三角网格非水密时失败。
另一种方法是使用广义绕数来确定符号。这是一个对于不干净的网格的鲁棒方法,但是效率较低,如果c的位置确定,则他的效率为
O
(
n
)
O(\sqrt n)
O(n)。
无论哪种情况,通过 igl::signed_distance 的接口都是:
// Choose type of signing to use
igl::SignedDistanceType type = SIGNED_DISTANCE_TYPE_PSEUDONORMAL;
igl::signed_distance(P,V,F,sign_type,S,I,C,N);
igl::point_mesh_squared_distance 的输出与上面相同,但现在 S 包含有符号(非平方)距离,并且额外的输出 N(仅在 type == SIGNED_DISTANCE_TYPE_PSEUDON 时设置)包含用于通过伪正态测试进行签名的法线。
示例704计算穿过兔子的切片上的带符号距离。
行进立方体
通常,3D 数据被捕获为空间上定义的标量场:
f
(
x
)
:
R
3
→
R
f(x):R^3\to R
f(x):R3→R。隐含在这个场中,标量场的等值面通常是显着的几何对象。值 v 处的等值面由
R
3
R^3
R3 中的所有点 x 组成,使得
f
(
x
)
=
v
f(x)=v
f(x)=v。几何处理的核心问题是将等值面提取为三角形网格,以进行进一步的基于网格的处理或可视化。 这称为等值线。
“行进立方体”42 是一种在规则格子(又名网格)上绘制等高线三线性函数 f 的著名方法。该方法的核心思想是根据单元每个顶点的函数值,使用从查找表中选择的预定义拓扑(也称为连通性)来绘制穿过每个单元的等值面(如果有的话)。该方法对网格中的所有单元(“立方体”)进行迭代(“行进”),并将最终的水密的网格缝合在一起。
在 libigl 中,igl::marching_cubes 根据在 nx × ny × nz 规则网格的顶点位置 GV 处采样的输入标量场 S 构造三角形网格 (V,F):
igl::marching_cubes(S,GV,nx,ny,nz,V,F);
(示例 705)对到输入网格的有符号距离进行采样(左),然后使用行进立方体重建表面以绘制 0 水平集的轮廓(中)。 为了进行比较,将这个带符号的距离场限制到指示函数和轮廓,其中显示了严重的混叠缺陷。
面片朝向
来自网络的模型有时会出现无方向性,因为每个三角形顶点的顺序不一致。 确定网格的一致面方向对于双面照明(例如,一侧有红色天鹅绒,另一侧有金丝的布料)和内外确定(例如,使用广义缠绕数)至关重要。
对于表示双面片的(开放)表面,libigl 提供了一个例程来强制网格的每个可定向补丁 (igl::orientable_patches) 内的方向一致:
igl::bfs_orient(F,FF,C);
这个简单的例程将在网格的每个面片上使用广度优先搜索,以在输出面 FF 中强制执行一致的面方向。
对于表示实体对象边界的(闭合或接近闭合)表面,libigl 提供了一个例程来重新定向面,以便顶点排序对应于顶点的逆时针排序,右手定则法线指向外侧。该方法 45[] 假设宇宙的大部分是空的。也就是说,空间中的大多数点都在固体物体的外部而不是内部。点是在曲面片上采样的。对于每个样本点,光线射入两个半球以计算每侧(距离加权)环境光遮挡的平均值。贴片的定向使得外侧被遮挡较少(较轻,即面向更多空隙空间)。
igl::embree::reorient_facets_raycast(V,F,FF,I);
布尔向量 I 揭示了 F 中的哪些行已在 FF 中翻转。
(示例 706)加载方向不一致的卡车模型(背面三角形显示为较暗)。 可定向贴片具有独特的颜色,然后面向外(左中)。 或者,每个单独的三角形都被视为一个“补丁”(右中)并独立向外定向。
扫描体积
移动固体物体 A 的扫掠体积 S 可以定义为空间中的任何点,使得在某个时刻该点位于固体内部。换句话说,它是刚体运动
f
(
t
)
f(t)
f(t)随时间变换的固体物体的并集:
S
=
⋃
t
∈
[
0
,
1
]
f
(
t
)
A
S=\bigcup_{t\in[0,1]}f(t)A
S=⋃t∈[0,1]f(t)A
由经历非正常旋转刚性运动三角形网格界定的实体扫掠体积,它的表面不是由三角形网格精确表示的表面:它将是分段规则表面。
要证实这一点,请考虑在执行螺旋运动时由单个边的线段扫过的表面。
这意味着,如果希望三角形网格的扫掠体积的表面经历刚性运动,并且希望输出是另一个三角形网格,那么将不得进行一定的近似误差。
考虑到这一点,计算近似扫掠体积的最简单方法是利用基于有符号距离的扫掠体积的替代定义:
S
=
{
p
∣
d
(
p
,
∂
S
)
<
0
}
=
{
p
∣
m
i
n
t
∈
[
0
,
1
]
d
(
p
,
f
(
t
)
∂
A
)
<
0
}
S=\{\bold p|d(\bold p,\partial S)< 0\}=\{\bold p|min_{t\in[0,1]}d(\bold p,f(t)\partial A)<0\}
S={p∣d(p,∂S)<0}={p∣mint∈[0,1]d(p,f(t)∂A)<0}
如果
∂
A
\partial A
∂A是一个三角网格,可以通过以下两种方式进行估计:1)以有限步长离散化时间
[
0
,
Δ
t
,
2
Δ
,
.
.
.
,
1
]
[0,\Delta t,2\Delta,...,1]
[0,Δt,2Δ,...,1]或是使用规则网格离散空间并使用网格值的三线性插值表示距离场。最后,使用移动立方体通过轮廓绘制来近似输出网格
∂
S
\partial S
∂S。
该方法类似于 Schroeder 等人在1994 年描述的方法。以及 Garg 等人在2006年与布尔运算结合使用的一个方法。
在 libigl 中,如果输入实体的表面由 (V,F) 表示,则调用以下代码后输出表面网格将为 (SV,SF):
igl::copyleft::swept_volume(V,F,num_time_steps,grid_size,isolevel,SV,SF);
等水平参数可以设置为零以近似精确的扫掠体积,大于零以近似扫掠体积的正偏移,或者小于零以近似负偏移。
(示例 707)计算经历刚性运动(金)的兔子模型的扫掠体积(银)的表面。
选取
使用鼠标拾取顶点和面在几何处理应用程序中非常常见。 虽然这看起来是一个简单的操作,但其实现并不简单。 Libigl 包含一个使用 Embree raycaster 解决此问题的函数。 示例 708 演示了其用法:
bool hit = igl::unproject_onto_mesh(
Vector2f(x,y),
F,
viewer.core.view * viewer.core.model,
viewer.core.proj,
viewer.core.viewport,
*ei,
fid,
bc);
此函数从视图平面沿视图方向投射光线。 变量x和y是鼠标屏幕坐标; view、model、proj分别是视图、模型和投影矩阵; viewport是OpenGL格式的视口; ei 包含由 Embree 构建的边界体积层次结构,fid 和 bc 分别是拾取的面和拾取位置的重心坐标。
(示例 708)通过射线投射进行拾取。 所选面的颜色为红色。
可扩展的局部内射映射
可扩展的局部内射映射算法允许在海量数据集上计算局部内射映射。 该算法与 ARAP 有许多相似之处,但使用重新加权方案来最小化任意失真能量,包括防止引入翻转的失真能量。
示例709包含三个演示:(1)大规模2D参数化的示例,(2)具有软约束的2D变形的示例,以及(3)具有软约束的3D变形的示例。 libigl 中的实现是独立的,并且依赖 Eigen 来求解全局步骤中使用的线性系统。 此处提供了依赖 Pardiso 的优化版本。
使用 SLIM 算法在 10 次迭代中计算具有 50k 个面的网格的局部单射参数化。
双射图的单纯复增强框架
单纯复数增强框架算法允许高效、鲁棒地计算双射映射。 该算法构建了一个脚手架结构,利用 SLIM 等高效的局部单射映射算法,保证了低失真的重叠自由映射,同时高效且可扩展。
示例710包含双射参数化骆驼网格的演示。
细分曲面
给定一个具有顶点 V 和面 F 的粗网格(又名笼),可以通过细分每个面来创建具有更多顶点和面的更高分辨率的网格。 也就是说,输入中的每个粗三角形都被许多较小的三角形替换。 Libigl 具有三种不同的方法来细分三角形网格。
“平面内”细分方法不会改变网格的点集或载体表面。 新顶点将添加到现有三角形的平面上,并且原始网格中幸存的顶点不会移动。
通过添加新面,细分算法会改变网格的组合。 组合数学的变化和定位高分辨率顶点的公式称为“细分规则”。
例如,在 igl::upsample 的平面内细分方法中,在每条边的中点处添加顶点:
v
a
b
=
1
2
(
v
a
+
v
b
)
v_{ab}=\frac{1}{2}(v_a+v_b)
vab=21(va+vb)并且每个三角面片
(
i
a
,
i
b
,
i
c
)
(i_a,i_b,i_c)
(ia,ib,ic)都会被四个新的三角面片
(
i
a
,
i
a
b
,
i
c
a
)
(i_a,i_{ab},i_{ca})
(ia,iab,ica),
(
i
b
,
i
b
c
,
i
a
b
)
(i_b,i_{bc},i_{ab})
(ib,ibc,iab),
(
i
a
b
,
i
b
c
,
i
c
a
)
(i_ab,i_{bc},i_{ca})
(iab,ibc,ica)以及
(
i
b
c
,
i
c
,
i
c
a
)
(i_{bc},i_c,i_{ca})
(ibc,ic,ica)该过程可以递归地应用,从而产生越来越细的网格。
igl::loop 的细分方法不在平面内。 细化网格的顶点被移动到其邻居的权重组合:网格在细化时被平滑。这种和其他平滑细分方法可以理解为样条曲线到曲面的推广。特别是,当我们考虑细分的递归应用的限制时,循环细分方法将收敛到
C
1
C^1
C1曲面。排除“不规则”或“异常”顶点(原笼的价数不等于6的顶点),曲面为
C
2
C^2
C2。igl::loop 和 igl::upsample 的组合(连通性和面数)是相同的:唯一的区别是顶点已在 igl::loop 中进行了平滑处理。
最后,libigl 还在 igl::false_barycentric_subdivision 中实现了平面内“假重心细分”的形式。该方法只是将每个三角形的重心
v
a
b
c
v_{abc}
vabc添加为新顶点,并将每个三角形替换为三个三角形
(
i
a
,
i
b
,
i
a
b
c
)
(i_a,i_{b},i_{abc})
(ia,ib,iabc),
(
i
b
,
i
c
,
i
a
b
c
)
(i_b,i_{c},i_{abc})
(ib,ic,iabc)以及
(
i
c
,
i
a
,
i
a
b
c
)
(i_c,i_{a},i_{abc})
(ic,ia,iabc)。与 igl::upsample 相比,此方法将创建内角越来越小的三角形,并且新顶点将以极端偏差对载体表面进行采样。
原始的粗网格和三种不同的细分方法:igl::upsample、igl::loop 和 igl::false_barycentric_subdivision。
数据平滑
定义在表面
Ω
\Omega
Ω 上的噪声函数
f
f
f可以使用能量最小化来平滑,该能量最小化平衡平滑项
E
S
E_S
ES 和二次拟合项:
u
=
a
r
g
m
i
n
u
α
E
s
(
u
)
+
(
1
−
α
)
∫
Ω
∣
∣
u
−
f
∣
∣
2
d
x
u=argmin_u\alpha E_s(u)+(1-\alpha)\int_{\Omega}||u-f||^2dx
u=argminuαEs(u)+(1−α)∫Ω∣∣u−f∣∣2dx
参数
α
\alpha
α 决定函数平滑的程度。
平滑能量的经典选择是具有零诺伊曼边界条件的函数的拉普拉斯能量,它是双调和能量的一种形式。它是使用余切拉普拉斯 L 和质量矩阵 M 构造的:
Q
L
=
L
′
∗
(
M
∖
L
)
QL = L'*(M\setminus L)
QL=L′∗(M∖L)构建的。然而,由于隐式零诺伊曼边界条件,如果 f 在边界处不具有零法向梯度,则函数行为在边界处会显着扭曲。
建议改用具有自然Hessian边界条件的双条和能量,它对应于矩阵
Q
H
=
H
′
∗
(
M
2
∖
H
)
QH=H'*(M2\setminus H)
QH=H′∗(M2∖H)的Hessian能量,其中 H 是有限元 Hessian 矩阵,M2 是堆叠质量矩阵 。 矩阵 H 和 QH 在 libigl 中分别实现为 igl::hessian 和 igl::hessian_energy。 示例 712 中给出了如何使用该函数的示例。
在下图中,可以清楚地看到具有零诺依曼边界条件的拉普拉斯能量与 Hessian 能量之间的差异:而第三张图像中的零诺依曼边界条件使函数的等值线偏置为垂直于边界,Hessian 能量 给出无偏结果。
(示例712)从左到右:甲虫网格上的函数、添加噪声的函数、使用拉普拉斯能量和零诺依曼边界条件进行平滑的结果以及使用Hessian能量进行平滑的结果。
形状预测
输入为为一组顶点
P
0
P_0
P0(不必是任何网格的一部分),以及一组约束
S
=
{
S
1
,
S
2
,
.
.
.
,
S
m
}
S=\{S_1,S_2,...,S_m\}
S={S1,S2,...,Sm},其中每个约束都是在
P
0
P_0
P0不同的稀疏子集上定义的。我们希望创建一组接近原始顶点集
P
0
P_0
P0的新点集
P
P
P(每个点都有对应的索引),同时遵守约束。可以实现其他目标,例如平滑度。 约束可以是非线性的,这使得问题非凸、困难并且没有保证的全局最优。 解决此类问题的一种非常流行的轻量级方法是局部全局迭代算法,包括以下两个步骤:
对于第k步的迭代:1.局部步骤:针对每个约束单独计算集合
P
κ
−
1
P_{\kappa-1}
Pκ−1到 S 的投影; 这意味着将多个约束中出现的每个点分解。 这可能是一个非线性操作,但如果约束稀疏,它就是一组许多小系统。2.全局步骤:将集合
P
κ
P_{\kappa}
Pκ整合到尽可能接近投影的碎片集合,并尽可能使用辅助目标函数。这会产生一个全局的二次目标函数。 此外,所得线性系统具有常数矩阵,因此可以进行预分解。
我们在 libigl 中实现的版本是Sofien Bouazi等描述的通用版本,它被包含在以下的两个头文件当中,<igl/shapeup.h> and <igl/shapeup_local_projections.h>。示例 713 中提供了实现规则性约束的演示(创建每个面尽可能规则的网格)。
局部步骤由 igl::shapeup_projection_function 类型的函数实例化。 全局步骤由两个函数完成:igl::shapeup_precomputation(),它预先计算全局步骤所需的矩阵;igl::shapeup_solve(),它根据初始解
P
0
P_0
P0和输入局部投影函数解决问题 。 数据结构 igl::ShapeUpData 包含运行算法所需的信息,并且可以配置; 例如,不言自明的变量 Maxiterations。
全局步骤最小化以下能量:
E
t
o
t
a
l
=
λ
s
h
a
p
e
E
s
h
a
p
e
+
λ
c
l
o
s
e
E
c
l
o
s
e
+
λ
s
m
o
o
t
h
E
s
m
o
o
t
h
E_{total}=\lambda_{shape}E_{shape}+\lambda_{close}E_{close}+\lambda_{smooth}E_{smooth}
Etotal=λshapeEshape+λcloseEclose+λsmoothEsmooth
其中系数
λ
\lambda
λ被编码在igl::ShapeUpData,并且可以在调用 igl::shapeup_precomputation() 之前更新。其中
E
s
h
a
p
e
E_{shape}
Eshape分量是积分能量(将
p
κ
p_{\kappa}
pκ你和到局部投影上)。
E
c
l
o
s
e
E_{close}
Eclose分量取决于位置约束,由参数b以及bc给出。而
E
s
m
o
o
t
h
E_{smooth}
Esmooth是一个可选的目标函数,用于最小化顶点之间的差异(在迪利克雷表达上的),通过参数E中的边进行编码。其中
E
c
l
o
s
e
E_{close}
Eclose以及
E
s
h
a
p
e
E_{shape}
Eshape分量都是各自通过wClose 以及wShape函数赋予权重的。
(示例 713)半隧道网格(左)已被优化为几乎完全规则(右)。色阶介于
[
0
,
0.05
]
[0,0.05]
[0,0.05],表征测量每个面角度与
9
0
∘
90^{\circ}
90∘的平均归一化偏差。
行进四面体
通常,3D 数据被捕获为空间上定义的标量场
f
(
x
)
:
R
3
→
R
f(x):R^3\to R
f(x):R3→R。隐含在这个场中,标量场的等值面通常是显着的几何对象。等值面在值v处由所有函数值
f
(
x
)
=
v
f(x)=v
f(x)=v的
R
3
R^3
R3中顶点x组成。几何处理的核心问题是将等值面提取为三角形网格,以进行进一步的基于网格的处理或可视化。 这称为等值线。
“行进四面体”是在 3D 单纯复形(也称为四面体网格)上等值线三线性函数 f 的著名方法。 该方法的核心思想是根据单元每个顶点的函数值,使用从查找表中选择的预定义拓扑(也称为连通性)来绘制穿过每个单元的等值面(如果有的话)。 该方法对复合体中的所有单元(“四面体”)进行迭代(“行进”),并将最终网格缝合在一起。
在 libigl 中,igl::marching_tets 构造一个三角形网格 (V,F),近似来自在四边形网格位置(TV、TT)顶点采样的输入标量场 S 的值 isovalue 的 iso-level 集:
igl::marching_tets(TV,TT,S, isovalue ,V,F);
隐式函数网格划分
待补充内容
快速测地距离近似的热方法
在上面的精确离散测地距离示例中,测地距离是精确计算的。这是一个十分耗时的操作:对于一个带有n条边的网格来说,它的计算复杂度为
O
(
n
2
l
o
g
(
n
)
)
O(n^2log(n))
O(n2log(n))。2013 年,Crane 等人提出了一种通过求解表面热方程、过滤结果然后通过求解泊松方程重建平滑解来更快地计算近似测地距离的方法。该方法从 Varadhan 的观察开始,两点 x 和 y 之间的测地距离 d 等于时间 t 后从 x 扩散到 y 的热量的对数的平方根。
d
(
x
,
y
)
=
l
i
m
t
→
0
−
4
t
l
o
g
k
t
,
x
(
y
)
d(x,y)=lim_{t\to 0}\sqrt{-4tlogk_{t,x}(y)}
d(x,y)=limt→0−4tlogkt,x(y)
其中
k
t
,
x
k_{t,x}
kt,x是热核心。我们可以将这个热扩散问题视为将热针放在 x 上,然后在 t 秒后测量点y的温度。
在三角形网格上,我们知道如何求解任意有限时间t内的热方程:
(
−
L
+
1
t
M
)
u
=
δ
x
(-\bold L +\frac{1}{t}\bold M)\bold u=\delta_x
(−L+t1M)u=δx
其中
L
∈
R
n
×
n
\bold L\in R^{n\times n}
L∈Rn×n是一个离散的拉普拉斯矩阵,
M
∈
R
n
×
n
\bold M\in R^{n\times n}
M∈Rn×n是一个离散的质量矩阵,
u
∈
R
n
\bold u\in R^{n}
u∈Rn是每个顶点的温度结果,而
δ
x
∈
R
n
\delta_x\in R^{n}
δx∈Rn是一个除了在源x处的顶点上为1而其他地方皆为零的一个向量。
如果我们有足够的数值准确度和精度,我们可以为一个很小的时间参数t简单地估计
−
4
t
l
o
g
u
\sqrt{-4tlog\ u}
−4tlog u。Crane 等人观察到的问题,是u的数值精度还远远不够。但是方向梯度
∇
u
\nabla \bold u
∇u却出奇的精确。因此,他们的想法是获取
u
\bold u
u的梯度,对这些向量进行归一化以获得梯度方向(单位向量)。 然后求解泊松方程,将这些方向积分为实际距离值。
该方法涉及到一个
n
×
n
n\times n
n×n矩阵的求逆运算(一个复杂度为
O
(
n
1....
)
O(n^{1....})
O(n1....)的操作),但如果使用 Cholesky 分解,则分解是预计算,即使测地距离的来源发生更改,也可以重复使用。 对于新的源,只需进行回代即可。
在 libigl 中,您可以使用此方法通过两个步骤计算从源顶点索引 gamma 列表到向量 D 的网格 (V,F) 的近似测地距离:
igl::HeatGeodesicsData<double> data;
igl::heat_geodesics_precompute(V,F,data);
...
igl::heat_geodesics_solve(data,gamma,D);
Example 716加载网格并计算从用户单击的位置开始的近似测地线距离。
隐式德劳内三角剖分
测地距离的原始热方法在规则的、无偏的网格上效果很好:即有限元余切和质量矩阵表现良好。 然而,对于质量较差的网格,此方法可能会显示任意较差的结果。 增加时间参数 t 可以减少这种不稳定性,但同时可以平滑所得的近似距离。
相反,一种改进途径是采用所谓的隐式 Delaunay 三角剖分离散拉普拉斯算子。
由于余切拉普拉斯算子仅取决于三角形网格的边长,因此该新算子将通过本质上翻转边并记录边长的变化来构造。 翻转边,直到每条边局部 Delaunay(即,其相应的余切权重为正)。
您可以使用以下方法计算 libigl 中网格 (V,F) 的隐式 Delaunay 三角剖分:
Eigen::MatrixXd l;
igl::edge_lengths(V,F,l);
Eigen::MatrixXd l_intrinsic;
Eigen::MatrixXi F_intrinsic;
igl::intrinsic_delaunay_triangulation(l,F,l_intrinsic,F_intrinsic);
请注意,未使用网格顶点位置 V,因为这是纯粹的内部操作。 该方法输入和输出为边长和三角形索引。
您可以直接使用以下方法构造内在 Delaunay 余切拉普拉斯矩阵:
Eigen::SparseMatrix<double> L;
igl::intrinsic_delaunay_cotmatrix(V,F,L);
最后,您可以通过以下方式使用此矩阵计算热测地线:
igl::HeatGeodesicsData<double> data;
data.use_intrinsic_delaunay = true;
igl::heat_geodesics_precompute(V,F,data);
...
igl::heat_geodesics_solve(data,gamma,D);
为三角形乱序集合以及云设计的快速绕组数
2018 年,Barill 等人演示了如何显着加快上述广义绕组数的计算。 三角形网格的广义缠绕数的原始定义也扩展到(定向)点云。
三角形乱序集合
对于三角形乱序集合,精确的分而治之方法具有理想的三角形数量对数
o
(
l
o
g
(
n
)
)
o(log(n))
o(log(n))级的复杂度。然而,该方法存在计算变为线性
o
(
n
)
o(n)
o(n)的故障模式。这可以通过类似于多体引力系统模拟或静电问题中使用的树方法来显着改善。结果是近似的,但更接近
o
(
l
o
g
(
n
)
o(log(n)
o(log(n)趋势并且具有更小的常数因子。
([例717]加载网格,对 1,000,000 个随机查询进行采样,然后丢弃给定模型之外的所有查询。)(images/bunny-fwn-soup.jpg)
计算三角形乱序集合的快速绕组数有两个步骤:构建树数据结构,然后在查询点进行评估。 在 libigl 中,编程如下:
igl::FastWindingNumberBVH fwn_bvh;
igl::fast_winding_number(V.cast<float>(),F,2,fwn_bvh);
Eigen::VectorXf W;
igl::fast_winding_number(fwn_bvh,2,Q.cast<float>(),W);
点云
对于点云,过程类似,但可能需要更多的预计算。 只要每个点都带有朝外的法线和面积值,就可以定义点云的缠绕数。 可以使用 libigl 函数来估计面积,该函数计算每个点 igl::copyleft::cgal::point_areas 的切线空间 Voronoi 图。 该函数又依赖于首先计算 k 个最近邻 igl::knn。 该函数和最终的绕数计算使用 libigl 的 igl::octree 作为包围体层次结构。 要估计具有法线 N 的点云 P 的面积,请使用:
// Build octree
std::vector<std::vector<int > > O_PI;
Eigen::MatrixXi O_CH;
Eigen::MatrixXd O_CN;
Eigen::VectorXd O_W;
igl::octree(P,O_PI,O_CH,O_CN,O_W);
Eigen::VectorXd A;
{
Eigen::MatrixXi I;
igl::knn(P,20,O_PI,O_CH,O_CN,O_W,I);
// CGAL is only used to help get point areas
igl::copyleft::cgal::point_areas(P,I,N,A);
}
然后可以计算查询列表 Q 的快速缠绕数:
Eigen::MatrixXd O_CM;
Eigen::VectorXd O_R;
Eigen::MatrixXd O_EC;
igl::fast_winding_number(P,N,A,O_PI,O_CH,2,O_CM,O_R,O_EC);
Eigen::VectorXd W;
igl::fast_winding_number(P,N,A,O_PI,O_CH,O_CM,O_R,O_EC,Q,2,W);