Libigl学习笔记——二、第七章——其他
- 第七章
- 7.1 网格统计 Mesh Statistics
- 7.2 广义绕组数 Generalized Winding Number
- 7.3 网格抽取 Mesh Decimation
- 7.4 有符号距离 Signed Distances
- 7.5 行进立方体 Marching Cubes
- 7.6 小平面方向 Facet Orientation
- 7.7 扫描卷 Swept Volume
- 7.8 拣选 Picking
- 7.9 可扩展的局部注入映射 Scalable Locally Injective Maps
- 7.10 双射映射的单纯复数增强框架 Simplicial Complex Augmentation Framework For Bijective Maps
- 7.11 细分曲面 Subdivision Surfaces
- 7.12 数据平滑 Data Smoothing
- 7.13 形状投影 Shapeup Projections
- 7.14 行进四面体 Marching Tetrahedra
- 7.15 隐式函数网格划分 Implicit Function Meshing
- 7.16 快速测地线距离近似的热法 Heat Method For Fast Geodesic Distance Approximation
- 7.17 汤和云的快速缠绕数 Fast Winding Number For Soups And Clouds
- 7.18 迭代最近点 Iterative Closest Point¶
- 7.19 分解图 Exploded View
- 7.20 蓝色噪声表面采样 Blue Noise Surface Sampling
- 7.21 矢量场平滑 Vector Field Smoothing
- 7.22 向量并行传输 Vector Parallel Transport
第七章
Libigl 包含各种各样的几何处理工具和函数，用于处理网格和与之相关的线性代数：本入门教程中讨论的内容太多了。我们在本章中列出了几个有趣的功能来重点介绍。
7.1 网格统计 Mesh Statistics
- Libigl 包含各种网格统计信息，包括面角度、面面积和奇异顶点的检测，奇异顶点是三角测量中邻域中或四边形中邻点多于或少于 4 个的顶点。
- 示例 Statistics 计算这些数量并进行基本的统计分析，以估计网格的等距和规则性：
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 的点。
- 第二行报告最小元素、最大元素和标准偏差的面积。这些数字按平均面积归一化，因此在上面的示例中，最大面积为 5.33 意味着最大的面比平均面大 5 倍。理想的各向同性网格的最小和最大面积都接近 1。
- 第三行测量面角，在完全规则的三角测量中，面角应接近 60 度（四边形为 90 度）。对于 FEM 而言，角度越接近 60 度，优化就越稳定。在这种情况下，很明显网格质量很差，如果用于求解偏微分方程，可能会导致伪影。
- 701案例完整展示：
#include <igl/doublearea.h>
#include <igl/internal_angles.h>
#include <igl/is_irregular_vertex.h>
#include <igl/readOBJ.h>
#include <igl/PI.h>
#include <Eigen/Core>
#include <iostream>
int main(int argc, char *argv[])
{
using namespace Eigen;
using namespace std;
MatrixXd V;
MatrixXi F;
igl::readOBJ(TUTORIAL_SHARED_PATH "/horse_quad.obj",V,F);
// Count the number of irregular vertices, the border is ignored
vector<bool> irregular = igl::is_irregular_vertex(V,F);
int vertex_count = V.rows();
int irregular_vertex_count =
std::count(irregular.begin(),irregular.end(),true);
double irregular_ratio = double(irregular_vertex_count)/vertex_count;
printf("Irregular vertices: \n%d/%d (%.2f%%)\n",
irregular_vertex_count,vertex_count, irregular_ratio*100);
// Compute areas, min, max and standard deviation
VectorXd area;
igl::doublearea(V,F,area);
area = area.array() / 2;
double area_avg = area.mean();
double area_min = area.minCoeff() / area_avg;
double area_max = area.maxCoeff() / area_avg;
double area_sigma = sqrt(((area.array()-area_avg)/area_avg).square().mean());
printf("Areas (Min/Max)/Avg_Area Sigma: \n%.2f/%.2f (%.2f)\n",
area_min,area_max,area_sigma);
// Compute per face angles, min, max and standard deviation
MatrixXd angles;
igl::internal_angles(V,F,angles);
angles = 360.0 * (angles/(2*igl::PI)); // Convert to degrees
double angle_avg = angles.mean();
double angle_min = angles.minCoeff();
double angle_max = angles.maxCoeff();
double angle_sigma = sqrt( (angles.array()-angle_avg).square().mean() );
printf("Angles in degrees (Min/Max) Sigma: \n%.2f/%.2f (%.2f)\n",
angle_min,angle_max,angle_sigma);
}
7.2 广义绕组数 Generalized Winding Number
- 封闭水密表面网格内部的四面体化问题是一个困难但适中的问题（参见我们的 Tetgen 包装器）。但是像TetGen这样的黑盒四边形网格器将拒绝来自多个连接组件的具有自相交、开放边界、非流形边的输入三角形网格。问题是双重的：自交存在矛盾的面约束，自交/开边界/非流形边使得在没有进一步假设的情况下确定内部和外部的问题变得不合适。
- 第一个问题可以通过“解决”所有自交叉点轻松解决。也就是说，对相交三角形进行网格划分，以便相交恰好发生在边和顶点处。这是使用 完成 igl::selfintersect 此操作的。
- TetGen通常可以将这个“已解决”网格的凸包四面体化，然后问题就变成了确定哪些t在输入网格内部，哪些在外部。也就是说，哪些应该保留，哪些应该删除。
- “广义绕组数”是一种确定麻烦网格内部和外部的可靠方法 40 。相对于某一 (V,F) 点 p∈R^3 的广义绕组数定义为标量函数：
- 其中 Ωfi ，在点 p 处由（第 i 个面） F 替换的 fi 立体角。这个立体角贡献是一个简单的闭式表达式，涉及 atan2 一些点积。
- 如果确实形成一个封闭的水密表面，那么 w§=1 如果位于内部， w§=0 如果 (V,F) p 位于外部 (V,F) (V,F) 。如果关闭但自身重叠，则 (V,F) w§ 是一个整数值，计算有多少（有符号）环绕 (V,F) p 。最后，如果不是闭合的或甚至不是流形（但至少是一致的定向），则 (V,F) w§ 平滑地趋向于1 p ，因为内部 (V,F) 更多，并且趋向于0 p ，因为外部更多。
- 702案例完整展示：
#include <igl/barycenter.h>
#include <igl/boundary_facets.h>
#include <igl/parula.h>
#include <igl/readMESH.h>
#include <igl/slice.h>
#include <igl/marching_tets.h>
#include <igl/winding_number.h>
#include <igl/opengl/glfw/Viewer.h>
#include <Eigen/Sparse>
#include <iostream>
Eigen::MatrixXd V,BC;
Eigen::VectorXd W;
Eigen::MatrixXi T,F,G;
double slice_z = 0.5;
enum OverLayType
{
OVERLAY_NONE = 0,
OVERLAY_INPUT = 1,
OVERLAY_OUTPUT = 2,
NUM_OVERLAY = 3,
} overlay = OVERLAY_NONE;
void update_visualization(igl::opengl::glfw::Viewer & viewer)
{
using namespace Eigen;
using namespace std;
Eigen::Vector4d plane(
0,0,1,-((1-slice_z)*V.col(2).minCoeff()+slice_z*V.col(2).maxCoeff()));
MatrixXd V_vis;
MatrixXi F_vis;
VectorXi J;
{
SparseMatrix<double> bary;
// Value of plane's implicit function at all vertices
const VectorXd IV =
(V.col(0)*plane(0) +
V.col(1)*plane(1) +
V.col(2)*plane(2)).array()
+ plane(3);
igl::marching_tets(V,T,IV,V_vis,F_vis,J,bary);
}
VectorXd W_vis;
igl::slice(W,J,W_vis);
MatrixXd C_vis;
// color without normalizing
igl::parula(W_vis,false,C_vis);
const auto & append_mesh = [&C_vis,&F_vis,&V_vis](
const Eigen::MatrixXd & V,
const Eigen::MatrixXi & F,
const RowVector3d & color)
{
F_vis.conservativeResize(F_vis.rows()+F.rows(),3);
F_vis.bottomRows(F.rows()) = F.array()+V_vis.rows();
V_vis.conservativeResize(V_vis.rows()+V.rows(),3);
V_vis.bottomRows(V.rows()) = V;
C_vis.conservativeResize(C_vis.rows()+F.rows(),3);
C_vis.bottomRows(F.rows()).rowwise() = color;
};
switch(overlay)
{
case OVERLAY_INPUT:
append_mesh(V,F,RowVector3d(1.,0.894,0.227));
break;
case OVERLAY_OUTPUT:
append_mesh(V,G,RowVector3d(0.8,0.8,0.8));
break;
default:
break;
}
viewer.data().clear();
viewer.data().set_mesh(V_vis,F_vis);
viewer.data().set_colors(C_vis);
viewer.data().set_face_based(true);
}
bool key_down(igl::opengl::glfw::Viewer& viewer, unsigned char key, int mod)
{
switch(key)
{
default:
return false;
case ' ':
overlay = (OverLayType)((1+(int)overlay)%NUM_OVERLAY);
break;
case '.':
slice_z = std::min(slice_z+0.01,0.99);
break;
case ',':
slice_z = std::max(slice_z-0.01,0.01);
break;
}
update_visualization(viewer);
return true;
}
int main(int argc, char *argv[])
{
using namespace Eigen;
using namespace std;
cout<<"Usage:"<<endl;
cout<<"[space] toggle showing input mesh, output mesh or slice "<<endl;
cout<<" through tet-mesh of convex hull."<<endl;
cout<<"'.'/',' push back/pull forward slicing plane."<<endl;
cout<<endl;
// Load mesh: (V,T) tet-mesh of convex hull, F contains facets of input
// surface mesh _after_ self-intersection resolution
igl::readMESH(TUTORIAL_SHARED_PATH "/big-sigcat.mesh",V,T,F);
// Compute barycenters of all tets
igl::barycenter(V,T,BC);
// Compute generalized winding number at all barycenters
cout<<"Computing winding number over all "<<T.rows()<<" tets..."<<endl;
igl::winding_number(V,F,BC,W);
// Extract interior tets
MatrixXi CT((W.array()>0.5).count(),4);
{
size_t k = 0;
for(size_t t = 0;t<T.rows();t++)
{
if(W(t)>0.5)
{
CT.row(k) = T.row(t);
k++;
}
}
}
// find bounary facets of interior tets
igl::boundary_facets(CT,G);
// boundary_facets seems to be reversed...
G = G.rowwise().reverse().eval();
// normalize
W = (W.array() - W.minCoeff())/(W.maxCoeff()-W.minCoeff());
// Plot the generated mesh
igl::opengl::glfw::Viewer viewer;
update_visualization(viewer);
viewer.callback_key_down = &key_down;
viewer.launch();
}
示例 702 计算具有孔和自相交点（金）的猫内部四面体网格的广义绕组数函数。银网是提取的内部Tet的表面，切片显示凸壳中所有Tet的缠绕数函数：蓝色（~ 0 ），绿色（~ 1），黄色（~2）
7.3 网格抽取 Mesh Decimation
- 网格简化或抽取的研究几乎与网格本身一样古老。给定一个包含太多三角形的高分辨率网格，请找到一个“近似”的低分辨率网格，其中包含更少的三角形。到目前为止，有各种不同的范式可以解决这个问题，最先进的方法相当先进。
- 一系列网格抽取方法通过从网格中连续删除元素来运行。特别是，霍普主张连续删除或塌陷边缘 39 []。该技术的通用形式是通过折叠单个边来构建从初始高分辨率网格到最低分辨率网格 M0 Mn 的 n 个网格序列：
- 霍普的原始方法和随后的后续工作提出了各种方法来选择下一个边缘在这个序列中崩溃。使用基于成本的范式，可以根据边缘的“成本”维护边缘的优先级队列（如果我删除此边缘，我的近似值会“差”多少？最便宜的边将折叠，相邻边的成本将更新。
- 为了保持拓扑（例如，如果网格组合为球体或圆环等），应该为坍缩会改变网格拓扑的边分配无限的成本。事实上，当且仅当坍缩边缘端点的共同邻居数量不正好是两个时，才会发生这种情况！
- 如果存在第三个共享顶点，则将删除另一个面，但将删除 2 条边。这可能会导致不必要的孔或非歧管“襟翼”。
有效的边折叠和无效的边折叠
- 还有一个一次性条件，即四面体的边缘不应折叠。
- 由于 libigl（有意地）不将其实现以动态网格数据结构（例如半边数据结构）为中心，因此对拓扑更改的支持有限。尽管如此，libigl 支持隔离边崩溃、边缘崩溃序列（每个在 O（log） 时间内）和基于优先级队列的抽取。
- 最简单的是 igl::decimation .通过调用
igl::decimate(V,F,1000,U,G);
- 网格将被抽取为新的网格 (V,F) (U,G) ，以便 G 具有大多数 1000 面。这将使用默认（朴素）条件来确定边折叠的成本和合并折点的放置。最短的边首先折叠，合并的顶点放置在边中点。
- 还可以提供函数句柄（ c++ lambda函数在这里很方便）， cost_and_placement 并 stopping_condition 分别用于确定边缘崩溃的成本/位置和停止条件。例如，上面的默认版本实现为：
igl::decimate(V,F,shortest_edge_and_midpoint,max_m,U,G);
- 其中 shortest_edge_and_midpoint ，将边的长度指定为成本，将其中点指定为合并的顶点放置，并 max_m 计算当前面数（有效折叠将计数减少 2）， true 如果计数低于 m=1000 。
- 还可以在抽取环路内更深入地刮擦并直接调用 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 个角对面的第 fth 面中（对于 o=0 边缘方向应该匹配，因为 o=1 方向相反）。
- 发生崩溃时， F 、 E 等矩阵的大小不会改变。相反，对应于“删除”面和边缘的行设置为特殊的常量值 IGL_COLLAPSE_EDGE_NULL 。这样做可以确保我们能够在真正恒定的时间 O（1） 内去除边缘。
- 方便. IGL_COLLAPSE_EDGE_NULL==0 这意味着大多数 OPENGL 风格的渲染 F 将简单地在第一个顶点绘制一堆 0 面积三角形。
- 以下内容将折叠第一条边并将其合并的折点放置在原点：
igl::collapse_edge(0,RowVector3d(0,0,0),V,F,E,EMAP,EF,EI);
- 如果有效，则相应地调整 、 V 、 F、 E、 EF 、 EI
- 这很强大，但水平很低。要围绕这一点构建一个抽取器，您需要跟踪哪些边要塌陷，哪些边接下来要塌陷。幸运的是，libigl 还公开了基于优先级队列的边缘折叠，其中包含用于调整成本和位置的功能句柄。
- 优先级队列被实现为（有序）集 Q 或（成本，边缘索引）对和迭代器列表，以便 Qit[e] 显示与eth边缘 Q 对应的迭代器 Qit 。展示位置存储在 #E 职位 C 列表中。调用以下内容时：
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演示了将这种基于优先级队列的方法与上面讨论的简单最短边中点成本/放置策略结合使用。
- 703示例完整展示：
#include <igl/circulation.h>
#include <igl/collapse_edge.h>
#include <igl/edge_flaps.h>
#include <igl/decimate.h>
#include <igl/shortest_edge_and_midpoint.h>
#include <igl/parallel_for.h>
#include <igl/read_triangle_mesh.h>
#include <igl/opengl/glfw/Viewer.h>
#include <Eigen/Core>
#include <iostream>
#include <set>
int main(int argc, char * argv[])
{
using namespace std;
using namespace Eigen;
using namespace igl;
cout<<"Usage: ./703_Decimation_bin [filename.(off|obj|ply)]"<<endl;
cout<<" [space] toggle animation."<<endl;
cout<<" 'r' reset."<<endl;
// Load a closed manifold mesh
string filename(TUTORIAL_SHARED_PATH "/fertility.off");
if(argc>=2)
{
filename = argv[1];
}
MatrixXd V,OV;
MatrixXi F,OF;
read_triangle_mesh(filename,OV,OF);
igl::opengl::glfw::Viewer viewer;
// Prepare array-based edge data structures and priority queue
VectorXi EMAP;
MatrixXi E,EF,EI;
igl::min_heap< std::tuple<double,int,int> > Q;
Eigen::VectorXi EQ;
// If an edge were collapsed, we'd collapse it to these points:
MatrixXd C;
int num_collapsed;
// Function to reset original mesh and data structures
const auto & reset = [&]()
{
F = OF;
V = OV;
edge_flaps(F,E,EMAP,EF,EI);
C.resize(E.rows(),V.cols());
VectorXd costs(E.rows());
// https://stackoverflow.com/questions/2852140/priority-queue-clear-method
// Q.clear();
Q = {};
EQ = Eigen::VectorXi::Zero(E.rows());
{
Eigen::VectorXd costs(E.rows());
igl::parallel_for(E.rows(),[&](const int e)
{
double cost = e;
RowVectorXd p(1,3);
shortest_edge_and_midpoint(e,V,F,E,EMAP,EF,EI,cost,p);
C.row(e) = p;
costs(e) = cost;
},10000);
for(int e = 0;e<E.rows();e++)
{
Q.emplace(costs(e),e,0);
}
}
num_collapsed = 0;
viewer.data().clear();
viewer.data().set_mesh(V,F);
viewer.data().set_face_based(true);
};
const auto &pre_draw = [&](igl::opengl::glfw::Viewer & viewer)->bool
{
// If animating then collapse 10% of edges
if(viewer.core().is_animating && !Q.empty())
{
bool something_collapsed = false;
// collapse edge
const int max_iter = std::ceil(0.01*Q.size());
for(int j = 0;j<max_iter;j++)
{
if(!collapse_edge(shortest_edge_and_midpoint,V,F,E,EMAP,EF,EI,Q,EQ,C))
{
break;
}
something_collapsed = true;
num_collapsed++;
}
if(something_collapsed)
{
viewer.data().clear();
viewer.data().set_mesh(V,F);
viewer.data().set_face_based(true);
}
}
return false;
};
const auto &key_down =
[&](igl::opengl::glfw::Viewer &viewer,unsigned char key,int mod)->bool
{
switch(key)
{
case ' ':
viewer.core().is_animating ^= 1;
break;
case 'R':
case 'r':
reset();
break;
default:
return false;
}
return true;
};
reset();
viewer.core().background_color.setConstant(1);
viewer.core().is_animating = true;
viewer.callback_key_down = key_down;
viewer.callback_pre_draw = pre_draw;
return viewer.launch();
}
示例703在生育模型上进行边缘塌陷
7.4 有符号距离 Signed Distances
- 在广义绕组数部分，我们研究了一种确定点是位于给定三角形汤网格内部还是外部的稳健方法。Libigl 通过加速有符号和无符号距离查询以及平面三角形网格和 3D 四面体网格的“元素内”查询来补充此算法。这些例程利用了 libigl 的通用轴对齐边界框层次结构 （ igl/AABB.h ）。此类是轻量级的，并且根据设计，不存储网格的副本（而是将其作为其成员函数的输入）。
7.4.1 点位置 Point Location
- 对于四面体网格，这对于“元素内”或“点位置”查询很有用：给定一个点 q∈R3 和一个四面体网格 (V,T) 确定四面体 q 所在的位置。这是在 libigl 中完成的，用于 tet 网格 V,T 和行中的查询点列表 igl::in_element() ， 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 可能属于多个四面体。在这些情况下，可以通过使用带有 a SparseMatrix 的 igl::in_element 重载作为输出来找到它们（而不仅仅是第一个）：
SparseMatrix<int> I;
igl::in_element(V,T,Q,tree,I);
- 现在每一 I 行都显示每个 tet 是否包含相应的行 Q ：表示 I(q,e)!=0 该点 q 在元素 e 中。
7.4.2 最近点 Closest Points
- 对于三角形网格，我们使用 AABB 树来加速点-网格最近点查询：给定网格 (V,F) 和查询点找到最近点 q∈R3 c∈(V,F) （其中 c 不一定是顶 (V,F) 点）。这是通过的三角形网格 V,F 和 via igl::point_mesh_squared_distance 行 P 中的点列表完成的：
VectorXd sqrD;
VectorXi I;
MatrixXd C;
igl::point_mesh_squared_distance(P,V,F,sqrD,I,C);
- 输出 sqrD 包含从每个点到 P 其最近点的（无符号）平方距离，该点 C 位于给定的元素 F 上 I （例如，可以使用从中恢复重心坐标 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);
7.4.3 符号距离 Signed Distance
- 最后，从最近的点或绕组编号开始，可以对这个距离进行签名。我们 igl::signed_distance 提供了两种签名方法：所谓的“伪正态检验” 37 []和广义绕组数 40 []。
- 伪正态检验（另见） igl::pseudonormal_test 假设输入网格是水密（闭合、非自相交、流形）网格。然后给定一个查询点和它的最近点 q c∈(V,F) ，它仔细选择一个外法线 n ， c 以便 sign(q−c)⋅n 揭示是在 q ： -1 内部 (V,F) 还是外部： +1。一旦找到，这是一个快速 O(1) 测试，但如果不是水密的， c 可能会 V,F 失败。
- 另一种方法是使用广义绕组编号来确定符号。这对于不干净的网格非常强大 V,F ，但速度较慢：类似于 O(√n) 一次 c 定位。
- 在任何一种情况下，接口通过 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 但现在包含有符号（非平方）距离，额外输出 N （仅在 时 type == SIGNED_DISTANCE_TYPE_PSEUDON 设置） S 包含用于使用伪正态测试签名的法线。
- 704案例完整展示：
#include <igl/cat.h>
#include <igl/edge_lengths.h>
#include <igl/parula.h>
#include <igl/per_edge_normals.h>
#include <igl/per_face_normals.h>
#include <igl/per_vertex_normals.h>
#include <igl/point_mesh_squared_distance.h>
#include <igl/readMESH.h>
#include <igl/signed_distance.h>
#include <igl/slice_mask.h>
#include <igl/marching_tets.h>
#include <igl/upsample.h>
#include <igl/opengl/glfw/Viewer.h>
#include <igl/writeOBJ.h>
#include <Eigen/Sparse>
#include <iostream>
Eigen::MatrixXd V;
Eigen::MatrixXi T,F;
igl::AABB<Eigen::MatrixXd,3> tree;
igl::FastWindingNumberBVH fwn_bvh;
Eigen::MatrixXd FN,VN,EN;
Eigen::MatrixXi E;
Eigen::VectorXi EMAP;
double max_distance = 1;
double slice_z = 0.5;
bool overlay = false;
bool useFastWindingNumber = false;
void update_visualization(igl::opengl::glfw::Viewer & viewer)
{
using namespace Eigen;
using namespace std;
Eigen::Vector4d plane(
0,0,1,-((1-slice_z)*V.col(2).minCoeff()+slice_z*V.col(2).maxCoeff()));
MatrixXd V_vis;
MatrixXi F_vis;
// Extract triangle mesh slice through volume mesh and subdivide nasty
// triangles
{
VectorXi J;
SparseMatrix<double> bary;
{
// Value of plane's implicit function at all vertices
const VectorXd IV =
(V.col(0)*plane(0) +
V.col(1)*plane(1) +
V.col(2)*plane(2)).array()
+ plane(3);
igl::marching_tets(V,T,IV,V_vis,F_vis,J,bary);
igl::writeOBJ("vis.obj",V_vis,F_vis);
}
while(true)
{
MatrixXd l;
igl::edge_lengths(V_vis,F_vis,l);
l /= (V_vis.colwise().maxCoeff() - V_vis.colwise().minCoeff()).norm();
const double max_l = 0.03;
if(l.maxCoeff()<max_l)
{
break;
}
Array<bool,Dynamic,1> bad = l.array().rowwise().maxCoeff() > max_l;
MatrixXi F_vis_bad, F_vis_good;
igl::slice_mask(F_vis,bad,1,F_vis_bad);
igl::slice_mask(F_vis,(bad!=true).eval(),1,F_vis_good);
igl::upsample(V_vis,F_vis_bad);
F_vis = igl::cat(1,F_vis_bad,F_vis_good);
}
}
// Compute signed distance
VectorXd S_vis;
if (!useFastWindingNumber)
{
VectorXi I;
MatrixXd N,C;
// Bunny is a watertight mesh so use pseudonormal for signing
signed_distance_pseudonormal(V_vis,V,F,tree,FN,VN,EN,EMAP,S_vis,I,C,N);
} else {
signed_distance_fast_winding_number(V_vis, V, F, tree, fwn_bvh, S_vis);
}
const auto & append_mesh = [&F_vis,&V_vis](
const Eigen::MatrixXd & V,
const Eigen::MatrixXi & F,
const RowVector3d & color)
{
F_vis.conservativeResize(F_vis.rows()+F.rows(),3);
F_vis.bottomRows(F.rows()) = F.array()+V_vis.rows();
V_vis.conservativeResize(V_vis.rows()+V.rows(),3);
V_vis.bottomRows(V.rows()) = V;
};
if(overlay)
{
append_mesh(V,F,RowVector3d(0.8,0.8,0.8));
}
viewer.data().clear();
viewer.data().set_mesh(V_vis,F_vis);
viewer.data().set_data(S_vis);
viewer.core().lighting_factor = overlay;
}
bool key_down(igl::opengl::glfw::Viewer& viewer, unsigned char key, int mod)
{
switch(key)
{
default:
return false;
case ' ':
overlay ^= true;
break;
case '.':
slice_z = std::min(slice_z+0.01,0.99);
break;
case ',':
slice_z = std::max(slice_z-0.01,0.01);
break;
case '1':
useFastWindingNumber = true;
break;
case '2':
useFastWindingNumber = false;
break;
}
update_visualization(viewer);
return true;
}
int main(int argc, char *argv[])
{
using namespace Eigen;
using namespace std;
cout<<"Usage:"<<endl;
cout<<"[space] toggle showing surface."<<endl;
cout<<"'.'/',' push back/pull forward slicing plane."<<endl;
cout<< "1/2 toggle between fast winding number (1) and pseudonormal (2) signing. \n";
cout<<endl;
// Load mesh: (V,T) tet-mesh of convex hull, F contains original surface
// triangles
igl::readMESH(TUTORIAL_SHARED_PATH "/bunny.mesh",V,T,F);
// Encapsulated call to point_mesh_squared_distance to determine bounds
{
VectorXd sqrD;
VectorXi I;
MatrixXd C;
igl::point_mesh_squared_distance(V,V,F,sqrD,I,C);
max_distance = sqrt(sqrD.maxCoeff());
}
// Fast winding and Pseudo normal depend on differnt AABB trees... We initialize both here.
// Pseudonormal setup...
// Precompute signed distance AABB tree
tree.init(V,F);
// Precompute vertex,edge and face normals
igl::per_face_normals(V,F,FN);
igl::per_vertex_normals(
V,F,igl::PER_VERTEX_NORMALS_WEIGHTING_TYPE_ANGLE,FN,VN);
igl::per_edge_normals(
V,F,igl::PER_EDGE_NORMALS_WEIGHTING_TYPE_UNIFORM,FN,EN,E,EMAP);
// fast winding number setup (just init fwn bvh)
igl::fast_winding_number(V, F, 2, fwn_bvh);
// Plot the generated mesh
igl::opengl::glfw::Viewer viewer;
update_visualization(viewer);
viewer.callback_key_down = &key_down;
viewer.data().show_lines = false;
viewer.launch();
}
示例 704 计算通过兔子的切片上的符号距离
7.5 行进立方体 Marching Cubes
- 通常，3D数据被捕获为在空间 f(x):R^3→R 上定义的标量场。潜伏在这个场中，标量场的等值面通常是突出的几何对象。值处的等值 v 面由 中的所有 R ^ 3 点 x 组成， f(x)=v 使得 .几何处理中的核心问题是将等值面提取为三角形网格，以便进一步进行基于网格的处理或可视化。这称为等轮廓。
- “行进立方体” 42 是在规则晶格（又名网格）上等轮廓三线性函数 f 的著名方法。此方法的核心思想是使用从查找表中选择的预定义拓扑（也称为连通性）勾勒出通过每个像元（如果有的话）的等值面轮廓，具体取决于像元每个顶点处的函数值。该方法遍历（“行进”）网格中的所有单元格（“立方体”），并将最终的水密网格拼接在一起。
- 在 libigl 中，从输入标量场构造一个三角形网格， igl::marching_cubes 该场 S 通过规则网格 (V,F) 在 nz a nx ny 的顶点位置 GV 采样：
igl::marching_cubes(S,GV,nx,ny,nz,V,F);
- 705案例完整展示：
#include <igl/marching_cubes.h>
#include <igl/signed_distance.h>
#include <igl/read_triangle_mesh.h>
#include <igl/voxel_grid.h>
#include <igl/opengl/glfw/Viewer.h>
#include <Eigen/Core>
#include <iostream>
int main(int argc, char * argv[])
{
using namespace Eigen;
using namespace std;
using namespace igl;
MatrixXi F;
MatrixXd V;
// Read in inputs as double precision floating point meshes
read_triangle_mesh(
TUTORIAL_SHARED_PATH "/armadillo.obj",V,F);
cout<<"Creating grid..."<<endl;
// number of vertices on the largest side
const int s = 100;
// create grid
MatrixXd GV;
Eigen::RowVector3i res;
igl::voxel_grid(V,0,s,1,GV,res);
// compute values
cout<<"Computing distances..."<<endl;
VectorXd S,B;
{
VectorXi I;
MatrixXd C,N;
signed_distance(GV,V,F,SIGNED_DISTANCE_TYPE_PSEUDONORMAL,S,I,C,N);
// Convert distances to binary inside-outside data --> aliasing artifacts
B = S;
for_each(B.data(),B.data()+B.size(),[](double& b){b=(b>0?1:(b<0?-1:0));});
}
cout<<"Marching cubes..."<<endl;
MatrixXd SV,BV;
MatrixXi SF,BF;
igl::marching_cubes(S,GV,res(0),res(1),res(2),0,SV,SF);
igl::marching_cubes(B,GV,res(0),res(1),res(2),0,BV,BF);
cout<<R"(Usage:
'1' Show original mesh.
'2' Show marching cubes contour of signed distance.
'3' Show marching cubes contour of indicator function.
)";
igl::opengl::glfw::Viewer viewer;
viewer.data().set_mesh(SV,SF);
viewer.callback_key_down =
[&](igl::opengl::glfw::Viewer & viewer, unsigned char key, int mod)->bool
{
switch(key)
{
default:
return false;
case '1':
viewer.data().clear();
viewer.data().set_mesh(V,F);
break;
case '2':
viewer.data().clear();
viewer.data().set_mesh(SV,SF);
break;
case '3':
viewer.data().clear();
viewer.data().set_mesh(BV,BF);
break;
}
viewer.data().set_face_based(true);
return true;
};
viewer.launch();
}
（示例705）对到输入网格（左）的符号距离进行采样，然后使用行进立方体重建表面以勾勒出0级集（中心）的轮廓。为了进行比较，将此有符号距离场固定到指示器函数并轮廓显示严重的混叠伪影。
7.6 小平面方向 Facet Orientation
- 来自 Web 的模型偶尔到达时没有定向，因为每个三角形顶点的顺序不一致。确定网格的一致刻面方向对于双面照明（例如，一侧为红色天鹅绒，另一侧为金丝的布）和内部和外部确定（例如，使用广义绕组编号）至关重要。
- 对于表示双面板的（开放）曲面，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案例完整展示：
#include <igl/read_triangle_mesh.h>
#include <igl/randperm.h>
#include <igl/orientable_patches.h>
#include <igl/slice.h>
#include <igl/hsv_to_rgb.h>
#include <igl/embree/reorient_facets_raycast.h>
#include <igl/opengl/glfw/Viewer.h>
#include <fstream>
#include <iostream>
#include <string>
#include <vector>
igl::opengl::glfw::Viewer viewer;
Eigen::MatrixXd V;
std::vector<Eigen::VectorXi> C(2);
std::vector<Eigen::MatrixXd> RGBcolors(2);
Eigen::MatrixXi F;
std::vector<Eigen::MatrixXi> FF(2);
bool is_showing_reoriented = false;
bool facetwise = false;
int main(int argc, char * argv[])
{
using namespace std;
cout<<R"(
Usage:
[space] Toggle between original and reoriented faces
F,f Toggle between patchwise and facetwise reorientation
S,s Scramble colors
)";
igl::read_triangle_mesh(TUTORIAL_SHARED_PATH "/truck.obj",V,F);
const auto & scramble_colors = []()
{
for(int pass = 0;pass<2;pass++)
{
Eigen::MatrixXi R;
igl::randperm(C[pass].maxCoeff()+1,R);
C[pass] = igl::slice(R,Eigen::MatrixXi(C[pass]));
Eigen::MatrixXd HSV(C[pass].rows(),3);
HSV.col(0) =
360.*C[pass].array().cast<double>()/(double)C[pass].maxCoeff();
HSV.rightCols(2).setConstant(1.0);
igl::hsv_to_rgb(HSV,RGBcolors[pass]);
}
viewer.data().set_colors(RGBcolors[facetwise]);
};
viewer.callback_key_pressed =
[&scramble_colors]
(igl::opengl::glfw::Viewer& /*viewer*/, unsigned int key, int mod)->bool
{
switch(key)
{
default:
return false;
case 'F':
case 'f':
{
facetwise = !facetwise;
break;
}
case 'S':
case 's':
{
scramble_colors();
return true;
}
case ' ':
{
is_showing_reoriented = !is_showing_reoriented;
break;
}
}
viewer.data().clear();
viewer.data().set_mesh(V,is_showing_reoriented?FF[facetwise]:F);
viewer.data().set_colors(RGBcolors[facetwise]);
return true;
};
// Compute patches
for(int pass = 0;pass<2;pass++)
{
Eigen::VectorXi I;
igl::embree::reorient_facets_raycast(
V,F,F.rows()*100,10,pass==1,false,false,I,C[pass]);
// apply reorientation
FF[pass].conservativeResize(F.rows(),F.cols());
for(int i = 0;i<I.rows();i++)
{
if(I(i))
{
FF[pass].row(i) = (F.row(i).reverse()).eval();
}else
{
FF[pass].row(i) = F.row(i);
}
}
}
viewer.data().set_mesh(V,is_showing_reoriented?FF[facetwise]:F);
viewer.data().set_face_based(true);
scramble_colors();
viewer.launch();
}
（示例 706）加载方向不一致的卡车模型（背面三角形显示较暗）。可定向的斑块具有独特的颜色，然后朝外（左中）。或者，每个单独的三角形都被视为一个“补丁”（中间右边）并独立地向外定向。
7.7 扫描卷 Swept Volume
-
移动固体物体 A 的扫描体积 S 可以定义为空间中的任何点，使得该点在某个时刻位于固体内部。换句话说，它是刚性运动 f(t) 随时间变化的固体物体的结合：
-
由三角形网格以非平凡旋转的刚性运动为界的实体的扫掠体积表面不是由三角形网格精确表示的表面：它将是一个分段规则表面。
-
要看到这一点，请考虑在执行螺钉运动时被单个边的线段扫描的曲面。
-
这意味着，如果我们想让三角形网格体的扫描体积表面经历刚性运动，并且我们希望输出是另一个三角形网格，那么我们必须对一定程度的近似误差感到满意。
-
考虑到这一点，计算近似扫描体积的最简单方法是利用基于符号距离的扫描体积的替代定义：
-
如果是一个三角形网格，那么 ∂A 我们可以通过以下方式近似：1）在有限的步长 [0,Δt,2Δt,…,1] 上离散时间和2）使用规则网格离散空间并使用网格值的三线性插值表示距离场。最后， ∂S 输出网格通过使用行进立方体 42 进行轮廓绘制来近似。
-
这种方法类似于Schroeder等人在1994年描述的方法 44 ，以及Garg等人在2016年与布尔运算结合使用的方法 38 。
-
在 libigl 中，如果输入实体的表面由 表示 (V,F) ，则输出表面网格将在调用后： (SV,SF)
igl::copyleft::swept_volume(V,F,num_time_steps,grid_size,isolevel,SV,SF);
- 该 isolevel 参数可以设置为零以近似确切的扫描体积，大于零表示近似扫描体积的正偏移，或小于零以近似于负偏移。
- 707案例完整展示：
#include <igl/read_triangle_mesh.h>
#include <igl/get_seconds.h>
#include <igl/material_colors.h>
#include <igl/swept_volume.h>
#include <igl/opengl/glfw/Viewer.h>
#include <igl/PI.h>
#include <Eigen/Core>
#include <iostream>
int main(int argc, char * argv[])
{
using namespace std;
using namespace igl;
Eigen::MatrixXi F,SF;
Eigen::MatrixXd V,SV,VT;
bool show_swept_volume = false;
// Define a rigid motion
const auto & transform = [](const double t)->Eigen::Affine3d
{
Eigen::Affine3d T = Eigen::Affine3d::Identity();
T.rotate(Eigen::AngleAxisd(t*2.*igl::PI,Eigen::Vector3d(0,1,0)));
T.translate(Eigen::Vector3d(0,0.125*cos(2.*igl::PI*t),0));
return T;
};
// Read in inputs as double precision floating point meshes
read_triangle_mesh(
TUTORIAL_SHARED_PATH "/bunny.off",V,F);
cout<<R"(Usage:
[space] Toggle between transforming original mesh and swept volume
)";
igl::opengl::glfw::Viewer viewer;
viewer.data().set_mesh(V,F);
viewer.data().set_face_based(true);
viewer.core().is_animating = !show_swept_volume;
const int grid_size = 50;
const int time_steps = 200;
const double isolevel = 0.1;
std::cerr<<"Computing swept volume...";
igl::swept_volume(
V,F,transform,time_steps,grid_size,isolevel,SV,SF);
std::cerr<<" finished."<<std::endl;
viewer.callback_pre_draw =
[&](igl::opengl::glfw::Viewer & viewer)->bool
{
if(!show_swept_volume)
{
Eigen::Affine3d T = transform(0.25*igl::get_seconds());
VT = V*T.matrix().block(0,0,3,3).transpose();
Eigen::RowVector3d trans = T.matrix().block(0,3,3,1).transpose();
VT = ( VT.rowwise() + trans).eval();
viewer.data().set_vertices(VT);
viewer.data().compute_normals();
}
return false;
};
viewer.callback_key_down =
[&](igl::opengl::glfw::Viewer & viewer, unsigned char key, int mod)->bool
{
switch(key)
{
default:
return false;
case ' ':
show_swept_volume = !show_swept_volume;
viewer.data().clear();
if(show_swept_volume)
{
viewer.data().set_mesh(SV,SF);
Eigen::Vector3d ambient = Eigen::Vector3d(SILVER_AMBIENT[0], SILVER_AMBIENT[1], SILVER_AMBIENT[2]);
Eigen::Vector3d diffuse = Eigen::Vector3d(SILVER_DIFFUSE[0], SILVER_DIFFUSE[1], SILVER_DIFFUSE[2]);
Eigen::Vector3d specular = Eigen::Vector3d(SILVER_SPECULAR[0], SILVER_SPECULAR[1], SILVER_SPECULAR[2]);
viewer.data().uniform_colors(ambient,diffuse,specular);
}
else
{
viewer.data().set_mesh(V,F);
}
viewer.core().is_animating = !show_swept_volume;
viewer.data().set_face_based(true);
break;
}
return true;
};
viewer.launch();
}
（示例707）计算兔子模型经历刚性运动（金）的扫描体积（银）的表面
7.8 拣选 Picking
- 使用鼠标拾取顶点和面在几何处理应用程序中非常常见。虽然这看起来很简单，但其实现并不简单。Libigl 包含一个使用 Embree 光线投射器解决此问题的函数。示例 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案例完整展示：
#include <igl/readOFF.h>
#include <igl/unproject_onto_mesh.h>
#include <igl/opengl/glfw/Viewer.h>
#include <iostream>
int main(int argc, char *argv[])
{
// Mesh with per-face color
Eigen::MatrixXd V, C;
Eigen::MatrixXi F;
// Load a mesh in OFF format
igl::readOFF(TUTORIAL_SHARED_PATH "/fertility.off", V, F);
// Initialize white
C = Eigen::MatrixXd::Constant(F.rows(),3,1);
igl::opengl::glfw::Viewer viewer;
viewer.callback_mouse_down =
[&V,&F,&C](igl::opengl::glfw::Viewer& viewer, int, int)->bool
{
int fid;
Eigen::Vector3f bc;
// Cast a ray in the view direction starting from the mouse position
double x = viewer.current_mouse_x;
double y = viewer.core().viewport(3) - viewer.current_mouse_y;
if(igl::unproject_onto_mesh(Eigen::Vector2f(x,y), viewer.core().view,
viewer.core().proj, viewer.core().viewport, V, F, fid, bc))
{
// paint hit red
C.row(fid)<<1,0,0;
viewer.data().set_colors(C);
return true;
}
return false;
};
std::cout<<R"(Usage:
[click] Pick face on shape
)";
// Show mesh
viewer.data().set_mesh(V, F);
viewer.data().set_colors(C);
viewer.data().show_lines = false;
viewer.launch();
}
（例708）通过光线投射进行采摘。所选面以红色着色
7.9 可扩展的局部注入映射 Scalable Locally Injective Maps
- 可扩展局部注入映射算法允许在大量数据集上计算本地注入映射 43 。该算法与ARAP有许多相似之处，但使用重新加权方案来最小化任意失真能量，包括那些防止引入翻转的能量。
- 示例 709 包含三个演示：（1） 大规模 2D 参数化示例，（2） 具有软约束的 2D 变形示例，以及 （3） 具有软约束的 3D 变形示例。libigl 中的实现是自包含的，并且依赖于特征来求解全局步骤中使用的线性系统。此处提供了依赖于Pardiso的优化版本。
- 709案例完整展示：
#include <iostream>
#include <igl/slim.h>
#include <igl/vertex_components.h>
#include <igl/readOBJ.h>
#include <igl/writeOBJ.h>
#include <igl/Timer.h>
#include <igl/boundary_loop.h>
#include <igl/map_vertices_to_circle.h>
#include <igl/harmonic.h>
#include <igl/MappingEnergyType.h>
#include <igl/serialize.h>
#include <igl/read_triangle_mesh.h>
#include <igl/opengl/glfw/Viewer.h>
#include <igl/flipped_triangles.h>
#include <igl/euler_characteristic.h>
#include <igl/barycenter.h>
#include <igl/adjacency_matrix.h>
#include <igl/is_edge_manifold.h>
#include <igl/doublearea.h>
#include <igl/cat.h>
#include <igl/PI.h>
#include <stdlib.h>
#include <string>
#include <vector>
using namespace std;
using namespace Eigen;
void check_mesh_for_issues(Eigen::MatrixXd& V, Eigen::MatrixXi& F);
void param_2d_demo_iter(igl::opengl::glfw::Viewer& viewer);
void get_soft_constraint_for_circle(Eigen::MatrixXd& V_o, Eigen::MatrixXi& F, Eigen::VectorXi& b, Eigen::MatrixXd& bc);
void soft_const_demo_iter(igl::opengl::glfw::Viewer& viewer);
void deform_3d_demo_iter(igl::opengl::glfw::Viewer& viewer);
void get_cube_corner_constraints(Eigen::MatrixXd& V_o, Eigen::MatrixXi& F, Eigen::VectorXi& b, Eigen::MatrixXd& bc);
void display_3d_mesh(igl::opengl::glfw::Viewer& viewer);
void int_set_to_eigen_vector(const std::set<int>& int_set, Eigen::VectorXi& vec);
Eigen::MatrixXd V;
Eigen::MatrixXi F;
bool first_iter = true;
igl::SLIMData sData;
igl::Timer timer;
double uv_scale_param;
enum DEMO_TYPE {
PARAM_2D,
SOFT_CONST,
DEFORM_3D
};
DEMO_TYPE demo_type;
bool key_down(igl::opengl::glfw::Viewer& viewer, unsigned char key, int modifier){
if (key == ' ') {
switch (demo_type) {
case PARAM_2D: {
param_2d_demo_iter(viewer);
break;
}
case SOFT_CONST: {
soft_const_demo_iter(viewer);
break;
}
case DEFORM_3D: {
deform_3d_demo_iter(viewer);
break;
}
default:
break;
}
}
return false;
}
void param_2d_demo_iter(igl::opengl::glfw::Viewer& viewer) {
if (first_iter) {
timer.start();
igl::read_triangle_mesh(TUTORIAL_SHARED_PATH "/face.obj", V, F);
check_mesh_for_issues(V,F);
cout << "\tMesh is valid!" << endl;
Eigen::MatrixXd uv_init;
Eigen::VectorXi bnd; Eigen::MatrixXd bnd_uv;
igl::boundary_loop(F,bnd);
igl::map_vertices_to_circle(V,bnd,bnd_uv);
igl::harmonic(V,F,bnd,bnd_uv,1,uv_init);
if (igl::flipped_triangles(uv_init,F).size() != 0) {
igl::harmonic(F,bnd,bnd_uv,1,uv_init); // use uniform laplacian
}
cout << "initialized parametrization" << endl;
sData.slim_energy = igl::MappingEnergyType::SYMMETRIC_DIRICHLET;
Eigen::VectorXi b; Eigen::MatrixXd bc;
slim_precompute(V,F,uv_init,sData, igl::MappingEnergyType::SYMMETRIC_DIRICHLET, b,bc,0);
uv_scale_param = 15 * (1./sqrt(sData.mesh_area));
viewer.data().set_mesh(V, F);
viewer.core().align_camera_center(V,F);
viewer.data().set_uv(sData.V_o*uv_scale_param);
viewer.data().compute_normals();
viewer.data().show_texture = true;
first_iter = false;
} else {
timer.start();
slim_solve(sData,1); // 1 iter
viewer.data().set_uv(sData.V_o*uv_scale_param);
}
cout << "time = " << timer.getElapsedTime() << endl;
cout << "energy = " << sData.energy << endl;
}
void soft_const_demo_iter(igl::opengl::glfw::Viewer& viewer) {
if (first_iter) {
igl::read_triangle_mesh(TUTORIAL_SHARED_PATH "/circle.obj", V, F);
check_mesh_for_issues(V,F);
cout << "\tMesh is valid!" << endl;
Eigen::MatrixXd V_0 = V.block(0,0,V.rows(),2);
Eigen::VectorXi b; Eigen::MatrixXd bc;
get_soft_constraint_for_circle(V_0,F,b,bc);
double soft_const_p = 1e5;
slim_precompute(V,F,V_0,sData,igl::MappingEnergyType::SYMMETRIC_DIRICHLET,b,bc,soft_const_p);
viewer.data().set_mesh(V, F);
viewer.core().align_camera_center(V,F);
viewer.data().compute_normals();
viewer.data().show_lines = true;
first_iter = false;
} else {
slim_solve(sData,1); // 1 iter
viewer.data().set_mesh(sData.V_o, F);
}
}
void deform_3d_demo_iter(igl::opengl::glfw::Viewer& viewer) {
if (first_iter) {
timer.start();
igl::readOBJ(TUTORIAL_SHARED_PATH "/cube_40k.obj", V, F);
Eigen::MatrixXd V_0 = V;
Eigen::VectorXi b; Eigen::MatrixXd bc;
get_cube_corner_constraints(V_0,F,b,bc);
double soft_const_p = 1e5;
sData.exp_factor = 5.0;
slim_precompute(V,F,V_0,sData,igl::MappingEnergyType::EXP_CONFORMAL,b,bc,soft_const_p);
//cout << "precomputed" << endl;
first_iter = false;
display_3d_mesh(viewer);
} else {
timer.start();
slim_solve(sData,1); // 1 iter
display_3d_mesh(viewer);
}
cout << "time = " << timer.getElapsedTime() << endl;
cout << "energy = " << sData.energy << endl;
}
void display_3d_mesh(igl::opengl::glfw::Viewer& viewer) {
MatrixXd V_temp; MatrixXi F_temp;
Eigen::MatrixXd Barycenters;
igl::barycenter(sData.V,sData.F,Barycenters);
//cout << "Barycenters.rows() = " << Barycenters.rows() << endl;
//double t = double((key - '1')+1) / 9.0;
double view_depth = 10.;
double t = view_depth/9.;
VectorXd v = Barycenters.col(2).array() - Barycenters.col(2).minCoeff();
v /= v.col(0).maxCoeff();
vector<int> s;
for (unsigned i=0; i<v.size();++i)
if (v(i) < t)
s.push_back(i);
V_temp.resize(s.size()*4,3);
F_temp.resize(s.size()*4,3);
for (unsigned i=0; i<s.size();++i){
V_temp.row(i*4+0) = sData.V_o.row(sData.F(s[i],0));
V_temp.row(i*4+1) = sData.V_o.row(sData.F(s[i],1));
V_temp.row(i*4+2) = sData.V_o.row(sData.F(s[i],2));
V_temp.row(i*4+3) = sData.V_o.row(sData.F(s[i],3));
F_temp.row(i*4+0) << (i*4)+0, (i*4)+1, (i*4)+3;
F_temp.row(i*4+1) << (i*4)+0, (i*4)+2, (i*4)+1;
F_temp.row(i*4+2) << (i*4)+3, (i*4)+2, (i*4)+0;
F_temp.row(i*4+3) << (i*4)+1, (i*4)+2, (i*4)+3;
}
viewer.data().set_mesh(V_temp,F_temp);
viewer.core().align_camera_center(V_temp,F_temp);
viewer.data().set_face_based(true);
viewer.data().show_lines = true;
}
int main(int argc, char *argv[]) {
cerr << "Press space for running an iteration." << std::endl;
cerr << "Syntax: " << argv[0] << " demo_number (1 to 3)" << std::endl;
cerr << "1. 2D unconstrained parametrization" << std::endl;
cerr << "2. 2D deformation with positional constraints" << std::endl;
cerr << "3. 3D mesh deformation with positional constraints" << std::endl;
demo_type = PARAM_2D;
if (argc == 2) {
switch (std::atoi(argv[1])) {
case 1: {
demo_type = PARAM_2D;
break;
} case 2: {
demo_type = SOFT_CONST;
break;
} case 3: {
demo_type = DEFORM_3D;
break;
}
default: {
cerr << "Wrong demo number - Please choose one between 1-3" << std:: endl;
exit(1);
}
}
}
// Launch the viewer
igl::opengl::glfw::Viewer viewer;
viewer.callback_key_down = &key_down;
// Disable wireframe
viewer.data().show_lines = false;
// Draw checkerboard texture
viewer.data().show_texture = false;
// First iteration
key_down(viewer, ' ', 0);
viewer.launch();
return 0;
}
void check_mesh_for_issues(Eigen::MatrixXd& V, Eigen::MatrixXi& F) {
Eigen::SparseMatrix<double> A;
igl::adjacency_matrix(F,A);
Eigen::MatrixXi C, Ci;
igl::vertex_components(A, C, Ci);
int connected_components = Ci.rows();
if (connected_components!=1) {
cout << "Error! Input has multiple connected components" << endl; exit(1);
}
int euler_char = igl::euler_characteristic(V, F);
if (euler_char!=1)
{
cout <<
"Error! Input does not have a disk topology, it's euler char is " <<
euler_char << endl;
exit(1);
}
bool is_edge_manifold = igl::is_edge_manifold(F);
if (!is_edge_manifold) {
cout << "Error! Input is not an edge manifold" << endl; exit(1);
}
Eigen::VectorXd areas; igl::doublearea(V,F,areas);
const double eps = 1e-14;
for (int i = 0; i < areas.rows(); i++) {
if (areas(i) < eps) {
cout << "Error! Input has zero area faces" << endl; exit(1);
}
}
}
void get_soft_constraint_for_circle(Eigen::MatrixXd& V_o, Eigen::MatrixXi& F, Eigen::VectorXi& b, Eigen::MatrixXd& bc) {
Eigen::VectorXi bnd;
igl::boundary_loop(F,bnd);
const int B_STEPS = 22; // constraint every B_STEPS vertices of the boundary
b.resize(bnd.rows()/B_STEPS);
bc.resize(b.rows(),2);
int c_idx = 0;
for (int i = B_STEPS; i < bnd.rows(); i+=B_STEPS) {
b(c_idx) = bnd(i);
c_idx++;
}
bc.row(0) = V_o.row(b(0)); // keep it there for now
bc.row(1) = V_o.row(b(2));
bc.row(2) = V_o.row(b(3));
bc.row(3) = V_o.row(b(4));
bc.row(4) = V_o.row(b(5));
bc.row(0) << V_o(b(0),0), 0;
bc.row(4) << V_o(b(4),0), 0;
bc.row(2) << V_o(b(2),0), 0.1;
bc.row(3) << V_o(b(3),0), 0.05;
bc.row(1) << V_o(b(1),0), -0.15;
bc.row(5) << V_o(b(5),0), +0.15;
}
void get_cube_corner_constraints(Eigen::MatrixXd& V, Eigen::MatrixXi& F, Eigen::VectorXi& b, Eigen::MatrixXd& bc) {
double min_x,max_x,min_y,max_y,min_z,max_z;
min_x = V.col(0).minCoeff(); max_x = V.col(0).maxCoeff();
min_y = V.col(1).minCoeff(); max_y = V.col(1).maxCoeff();
min_z = V.col(2).minCoeff(); max_z = V.col(2).maxCoeff();
// get all cube corners
b.resize(8,1); bc.resize(8,3);
int x;
for (int i = 0; i < V.rows(); i++) {
if (V.row(i) == Eigen::RowVector3d(min_x,min_y,min_z)) b(0) = i;
if (V.row(i) == Eigen::RowVector3d(min_x,min_y,max_z)) b(1) = i;
if (V.row(i) == Eigen::RowVector3d(min_x,max_y,min_z)) b(2) = i;
if (V.row(i) == Eigen::RowVector3d(min_x,max_y,max_z)) b(3) = i;
if (V.row(i) == Eigen::RowVector3d(max_x,min_y,min_z)) b(4) = i;
if (V.row(i) == Eigen::RowVector3d(max_x,max_y,min_z)) b(5) = i;
if (V.row(i) == Eigen::RowVector3d(max_x,min_y,max_z)) b(6) = i;
if (V.row(i) == Eigen::RowVector3d(max_x,max_y,max_z)) b(7) = i;
}
// get all cube edges
std::set<int> cube_edge1; Eigen::VectorXi cube_edge1_vec;
for (int i = 0; i < V.rows(); i++) {
if ((V(i,0) == min_x && V(i,1) == min_y)) {
cube_edge1.insert(i);
}
}
Eigen::VectorXi edge1;
int_set_to_eigen_vector(cube_edge1, edge1);
std::set<int> cube_edge2; Eigen::VectorXi edge2;
for (int i = 0; i < V.rows(); i++) {
if ((V(i,0) == max_x && V(i,1) == max_y)) {
cube_edge2.insert(i);
}
}
int_set_to_eigen_vector(cube_edge2, edge2);
b = igl::cat(1,edge1,edge2);
std::set<int> cube_edge3; Eigen::VectorXi edge3;
for (int i = 0; i < V.rows(); i++) {
if ((V(i,0) == max_x && V(i,1) == min_y)) {
cube_edge3.insert(i);
}
}
int_set_to_eigen_vector(cube_edge3, edge3);
b = igl::cat(1,b,edge3);
std::set<int> cube_edge4; Eigen::VectorXi edge4;
for (int i = 0; i < V.rows(); i++) {
if ((V(i,0) == min_x && V(i,1) == max_y)) {
cube_edge4.insert(i);
}
}
int_set_to_eigen_vector(cube_edge4, edge4);
b = igl::cat(1,b,edge4);
bc.resize(b.rows(),3);
Eigen::Matrix3d m; m = Eigen::AngleAxisd(0.3 * igl::PI, Eigen::Vector3d(1./sqrt(2.),1./sqrt(2.),0.)/*Eigen::Vector3d::UnitX()*/);
int i = 0;
for (; i < cube_edge1.size(); i++) {
Eigen::RowVector3d edge_rot_center(min_x,min_y,(min_z+max_z)/2.);
bc.row(i) = (V.row(b(i)) - edge_rot_center) * m + edge_rot_center;
}
for (; i < cube_edge1.size() + cube_edge2.size(); i++) {
Eigen::RowVector3d edge_rot_center(max_x,max_y,(min_z+max_z)/2.);
bc.row(i) = (V.row(b(i)) - edge_rot_center) * m.transpose() + edge_rot_center;
}
for (; i < cube_edge1.size() + cube_edge2.size() + cube_edge3.size(); i++) {
bc.row(i) = 0.75*V.row(b(i));
}
for (; i < b.rows(); i++) {
bc.row(i) = 0.75*V.row(b(i));
}
}
void int_set_to_eigen_vector(const std::set<int>& int_set, Eigen::VectorXi& vec) {
vec.resize(int_set.size()); int idx = 0;
for(auto f : int_set) {
vec(idx) = f; idx++;
}
}
使用 SLIM 算法在 10 次迭代中计算具有 50k 面的网格的局部注入参数化。
7.10 双射映射的单纯复数增强框架 Simplicial Complex Augmentation Framework For Bijective Maps
- 单纯复数增强框架 49 算法允许高效、鲁棒地计算双射映射。该算法构建了一个脚手架结构，以利用SLIM等高效的局部注入映射算法，保证了低失真的重叠自由映射，同时具有高效和可扩展性。
- 示例 710 包含对驼峰网格进行双射参数化的演示。
- 710示例完整展示：
......
示例 710 包含对驼峰网格进行双射参数化的演示
7.11 细分曲面 Subdivision Surfaces
-
给定一个带有顶点和面的粗网格（又名笼子），可以通过细分每个面来创建具有更多顶点 V 和面 F 的更高分辨率网格。也就是说，输入中的每个粗三角形都被许多较小的三角形所取代。Libigl 有三种不同的方法来细分三角形网格。
-
“平面内”细分方法不会改变网格的点集或载体表面。在现有三角形的平面上添加新顶点，并且不会移动从原始网格中幸存下来的顶点。
-
通过添加新面，细分算法会改变网格的组合。组合数学的变化和高分辨率顶点定位公式称为“细分规则”。
-
-
-
-
711 案例完整展示：
#include <igl/read_triangle_mesh.h>
#include <igl/loop.h>
#include <igl/upsample.h>
#include <igl/false_barycentric_subdivision.h>
#include <igl/opengl/glfw/Viewer.h>
#include <Eigen/Core>
#include <iostream>
int main(int argc, char * argv[])
{
using namespace std;
using namespace igl;
Eigen::MatrixXi OF,F;
Eigen::MatrixXd OV,V;
bool show_swept_volume = false;
read_triangle_mesh(
TUTORIAL_SHARED_PATH "/decimated-knight.off",OV,OF);
V = OV;
F = OF;
cout<<R"(Usage:
1 Restore Original mesh
2 Apply In-plane upsampled mesh
3 Apply Loop subdivided mesh
4 Apply False barycentric subdivision
)";
igl::opengl::glfw::Viewer viewer;
viewer.data().set_mesh(V,F);
viewer.data().set_face_based(true);
viewer.callback_key_down =
[&](igl::opengl::glfw::Viewer & viewer, unsigned char key, int mod)->bool
{
switch(key)
{
default:
return false;
case '1':
{
V = OV;
F = OF;
break;
}
case '2':
{
igl::upsample( Eigen::MatrixXd(V), Eigen::MatrixXi(F), V,F);
break;
}
case '3':
{
igl::loop( Eigen::MatrixXd(V), Eigen::MatrixXi(F), V,F);
break;
}
case '4':
{
igl::false_barycentric_subdivision(
Eigen::MatrixXd(V),Eigen::MatrixXi(F),V,F);
break;
}
}
viewer.data().clear();
viewer.data().set_mesh(V,F);
viewer.data().set_face_based(true);
return true;
};
viewer.launch();
}
原始粗网格和三种不同的细分方法： igl::upsample 、 igl::loop 和 igl::false_barycentric_subdivision
7.12 数据平滑 Data Smoothing
- 曲面上 Ω 定义的噪声函数 f 可以使用能量最小化来平滑，该能量最小化平衡平滑项与二次拟合项 ES ：
- 该参数 α 确定函数平滑的力度。
- 平滑度能量的经典选择是诺依曼边界条件为零的函数的拉普拉斯能量，它是双谐能量的一种形式。它是使用余切拉普拉斯量 L 和质量矩阵 M 构造的 QL = L’*(M\L) ：然而，由于隐含的零诺依曼边界条件，如果在边界处没有零法向梯度，则 f 函数行为在边界处明显扭曲。
- 其中建议使用具有自然黑森边界条件的双调和能量，这与具有矩阵的黑森能量相对应，其中 51 H 是有限元黑森， M2 是堆叠质量矩阵 QH = H’*(M2\H) 。矩阵 H 和 QH 分别以 libigl igl::hessian igl::hessian_energy 实现。示例 712 中给出了如何使用该函数的示例。
- 在下图中，可以清楚地看到具有零诺依曼边界条件的拉普拉斯能量与黑森能量之间的差异：第三幅图像中的零诺依曼边界条件使函数的等值线垂直于边界，而黑森能量给出了无偏的结果。
- 712案例完整展示：
#include <igl/read_triangle_mesh.h>
#include <igl/hessian_energy.h>
#include <igl/curved_hessian_energy.h>
#include <igl/massmatrix.h>
#include <igl/cotmatrix.h>
#include <igl/isolines_map.h>
#include <igl/parula.h>
#include <igl/vertex_components.h>
#include <igl/remove_unreferenced.h>
#include <igl/opengl/glfw/Viewer.h>
#include <igl/heat_geodesics.h>
#include <Eigen/Core>
#include <Eigen/SparseCholesky>
#include <iostream>
#include <set>
#include <limits>
#include <stdlib.h>
int main(int argc, char * argv[])
{
typedef Eigen::SparseMatrix<double> SparseMat;
srand(57);
//Read our mesh
Eigen::MatrixXd V;
Eigen::MatrixXi F;
if(!igl::read_triangle_mesh(
argc>1?argv[1]: TUTORIAL_SHARED_PATH "/beetle.off",V,F)) {
std::cout << "Failed to load mesh." << std::endl;
}
//Constructing an exact function to smooth
igl::HeatGeodesicsData<double> hgData;
igl::heat_geodesics_precompute(V, F, hgData);
Eigen::VectorXd heatDist;
Eigen::VectorXi gamma(1); gamma << 1947; //1631;
igl::heat_geodesics_solve(hgData, gamma, heatDist);
Eigen::VectorXd zexact =
0.1*(heatDist.array() + (-heatDist.maxCoeff())).pow(2)
+ 3*V.block(0,1,V.rows(),1).array().cos();
//Make the exact function noisy
const double s = 0.1*(zexact.maxCoeff() - zexact.minCoeff());
Eigen::VectorXd znoisy = zexact + s*Eigen::VectorXd::Random(zexact.size());
//Constructing the squared Laplacian and squared Hessian energy
SparseMat L, M;
igl::cotmatrix(V, F, L);
igl::massmatrix(V, F, igl::MASSMATRIX_TYPE_BARYCENTRIC, M);
Eigen::SimplicialLDLT<SparseMat> solver(M);
SparseMat MinvL = solver.solve(L);
SparseMat QL = L.transpose()*MinvL;
SparseMat QH;
igl::hessian_energy(V, F, QH);
SparseMat QcH;
igl::curved_hessian_energy(V, F, QcH);
//Solve to find Laplacian-smoothed Hessian-smoothed, and
// curved-Hessian-smoothed solutions
const double al = 3e-7;
Eigen::SimplicialLDLT<SparseMat> lapSolver(al*QL + (1.-al)*M);
Eigen::VectorXd zl = lapSolver.solve(al*M*znoisy);
const double ah = 2e-7;
Eigen::SimplicialLDLT<SparseMat> hessSolver(ah*QH + (1.-ah)*M);
Eigen::VectorXd zh = hessSolver.solve(ah*M*znoisy);
const double ach = 3e-7;
Eigen::SimplicialLDLT<SparseMat> curvedHessSolver(al*QcH + (1.-ach)*M);
Eigen::VectorXd zch = curvedHessSolver.solve(ach*M*znoisy);
//Viewer that shows all functions: zexact, znoisy, zl, zh
igl::opengl::glfw::Viewer viewer;
viewer.data().set_mesh(V,F);
viewer.data().show_lines = false;
viewer.callback_key_down =
[&](igl::opengl::glfw::Viewer & viewer, unsigned char key, int mod)->bool
{
//Graduate result to show isolines, then compute color matrix
const Eigen::VectorXd* z;
switch(key) {
case '1':
z = &zexact;
break;
case '2':
z = &znoisy;
break;
case '3':
z = &zl;
break;
case '4':
z = &zh;
break;
case '5':
z = &zch;
break;
default:
return false;
}
viewer.data().set_data(*z);
return true;
};
std::cout << R"(Smoothing a noisy function.
Usage:
1 Show original function
2 Show noisy function
3 Biharmonic smoothing (zero Neumann boundary)
4 Biharmonic smoothing (natural planar Hessian boundary)
5 Biharmonic smoothing (natural curved Hessian boundary)
)";
Eigen::MatrixXd CM;
igl::parula(Eigen::VectorXd::LinSpaced(21,0,1).eval(),false,CM);
igl::isolines_map(Eigen::MatrixXd(CM),CM);
viewer.data().set_colormap(CM);
viewer.data().set_data(znoisy);
viewer.launch();
//Constructing a step function to smooth
Eigen::VectorXd zstep = Eigen::VectorXd::Zero(V.rows());
for(int i=0; i<V.rows(); ++i) {
zstep(i) = V(i,2)<-0.25 ? 1. : (V(i,2)>0.31 ? 2. : 0);
}
//Smooth that function
const double sl = 2e-5;
Eigen::SimplicialLDLT<SparseMat> stepLapSolver(sl*QL + (1.-sl)*M);
Eigen::VectorXd stepzl = stepLapSolver.solve(al*M*zstep);
const double sh = 6e-6;
Eigen::SimplicialLDLT<SparseMat> stepHessSolver(sh*QH + (1.-sh)*M);
Eigen::VectorXd stepzh = stepHessSolver.solve(ah*M*zstep);
const double sch = 2e-5;
Eigen::SimplicialLDLT<SparseMat> stepCurvedHessSolver(sl*QcH + (1.-sch)*M);
Eigen::VectorXd stepzch = stepCurvedHessSolver.solve(ach*M*zstep);
//Display functions
igl::opengl::glfw::Viewer viewer2;
viewer2.data().set_mesh(V,F);
viewer2.data().show_lines = false;
viewer2.callback_key_down =
[&](igl::opengl::glfw::Viewer & viewer, unsigned char key, int mod)->bool
{
//Graduate result to show isolines, then compute color matrix
const Eigen::VectorXd* z;
switch(key) {
case '1':
z = &zstep;
break;
case '2':
z = &stepzl;
break;
case '3':
z = &stepzh;
break;
case '4':
z = &stepzch;
break;
default:
return false;
}
viewer.data().set_data(*z);
return true;
};
std::cout << R"(Smoothing a step function.
Usage:
1 Show step function
2 Biharmonic smoothing (zero Neumann boundary)
3 Biharmonic smoothing (natural planar Hessian boundary)
4 Biharmonic smoothing (natural curved Hessian boundary)
)";
viewer2.data().set_colormap(CM);
viewer2.data().set_data(zstep);
viewer2.launch();
return 0;
}
例 712）从左到右：甲虫网格上的函数，增加噪声的函数，拉普拉斯能和零诺依曼边界条件平滑的结果，以及黑森能量平滑的结果
7.13 形状投影 Shapeup Projections
- 我们的输入是一组点 P0 （不一定是任何网格的一部分）和一组约束，其中每个约束 S={S1,S2,…Sm} 都定义在 的不同 P0 且稀疏的子集上。我们希望创建一组接近原始集合 P0 的新点（每个点 P 都有相应的索引），同时遵守约束。可以采用其他目标，例如平滑度。约束可以是非线性的，这使得问题非凸、困难且没有保证的全局最优。解决此类问题的一种非常流行的轻量级方法是局部-全局迭代算法，包括以下两个步骤：
- 对于迭代：1.局部步骤 k ：计算集合 Pk−1 的投影到 S ，每个约束单独计算;这意味着将出现在多个约束中的每个点分开。这可以是一个非线性操作，但如果约束很稀疏，它是一组许多小系统。2.全局步骤：将集合 Pk 积分为尽可能接近投影的碎片集，并可能具有辅助目标函数。这就产生了一个全局但二次的目标函数。此外，得到的线性系统具有常数矩阵，因此可以进行预分解。
- 我们在 libigl 中实现的版本是 描述 27 的一般版本，位于两个文件中： <igl/shapeup.h> 和<igl/shapeup_local_projections.h> 。实现规则性约束（创建每个面尽可能规则的网格）的演示如例 713 所示。
- 本地步骤由 类型的 igl::shapeup_projection_function 函数实例化。全局步长由两个函数完成： igl::shapeup_precomputation() ，它根据初始解 P0 和输入局部投影函数预先计算全局步长所需的矩阵，以及igl::shapeup_solve() ，它解决问题。数据结构 igl::ShapeUpData 包含运行算法所需的信息，并且可以配置;例如，不言自明的变量 Maxiterations 。
- 全局步骤可最大程度地减少以下能量：
- 其中系数 λ 以 编码 igl::ShapeUpData ，并且可以在调用 igl::shapeup_precomputation() 之前更新。该 Eshape 分量是积分能量（拟合 Pk 局部投影）。该 Eclose 组件遵守位置约束，由 b 和 bc 参数给出。该 Esmooth 组件是一个可选的目标函数，为了最小化点之间的差异（在狄利克雷意义上），通过参数 E 中的“边”进行编码。两者 Eshape 也 Eclose 分别由 和 wShape 函数参数加 wClose 权。
- 713案例完整展示：
#include <igl/avg_edge_length.h>
#include <igl/barycenter.h>
#include <igl/jet.h>
#include <igl/shapeup.h>
#include <igl/quad_planarity.h>
#include <igl/readDMAT.h>
#include <igl/readOFF.h>
#include <igl/slice.h>
#include <igl/opengl/glfw/Viewer.h>
#include <igl/PI.h>
#include <vector>
#include <cstdlib>
// Quad mesh loaded
Eigen::MatrixXd VQC;
Eigen::MatrixXi FQC;
Eigen::MatrixXi E;
Eigen::MatrixXi FQCtri;
Eigen::MatrixXd PQC0, PQC1, PQC2, PQC3;
// Euclidean-regular quad mesh
Eigen::MatrixXd VQCregular;
Eigen::MatrixXi FQCtriregular;
Eigen::MatrixXd PQC0regular, PQC1regular, PQC2regular, PQC3regular;
igl::ShapeupData su_data;
// Scale for visualizing the fields
double global_scale; //TODO: not used
void quadAngleRegularity(const Eigen::MatrixXd& V, const Eigen::MatrixXi& Q, Eigen::VectorXd& angleRegularity)
{
angleRegularity.conservativeResize(Q.rows());
angleRegularity.setZero();
for (int i=0;i<Q.rows();i++){
for (int j=0;j<4;j++){
Eigen::RowVectorXd v21=(V.row(Q(i,j))-V.row(Q(i,(j+1)%4))).normalized();
Eigen::RowVectorXd v23=(V.row(Q(i,(j+2)%4))-V.row(Q(i,(j+1)%4))).normalized();
angleRegularity(i)+=(abs(acos(v21.dot(v23))-igl::PI/2.0)/(igl::PI/2.0))/4.0;
}
}
}
bool key_down(igl::opengl::glfw::Viewer& viewer, unsigned char key, int modifier)
{
using namespace std;
using namespace Eigen;
// Plot the original quad mesh
if (key == '1')
{
viewer.data().clear();
// Draw the triangulated quad mesh
viewer.data().set_mesh(VQC, FQCtri);
// Assign a color to each quad that corresponds to the average deviation of each angle from pi/2
VectorXd angleRegularity(FQC.rows());
quadAngleRegularity( VQC, FQC, angleRegularity);
MatrixXd Ct;
igl::jet(angleRegularity, 0.0, 0.05, Ct);
MatrixXd C(FQCtri.rows(),3);
C << Ct, Ct;
viewer.data().set_colors(C);
// Plot a line for each edge of the quad mesh
viewer.data().add_edges(PQC0, PQC1, Eigen::RowVector3d(0,0,0));
viewer.data().add_edges(PQC1, PQC2, Eigen::RowVector3d(0,0,0));
viewer.data().add_edges(PQC2, PQC3, Eigen::RowVector3d(0,0,0));
viewer.data().add_edges(PQC3, PQC0, Eigen::RowVector3d(0,0,0));
}
// Plot the planarized quad mesh
if (key == '2')
{
viewer.data().clear();
// Draw the triangulated quad mesh
viewer.data().set_mesh(VQCregular, FQCtri);
// Assign a color to each quad that corresponds to its planarity
VectorXd angleRegularity(FQC.rows());
quadAngleRegularity( VQCregular, FQC, angleRegularity);
MatrixXd Ct;
igl::jet(angleRegularity, 0, 0.05, Ct);
MatrixXd C(FQCtri.rows(),3);
C << Ct, Ct;
viewer.data().set_colors(C);
// Plot a line for each edge of the quad mesh
viewer.data().add_edges(PQC0regular, PQC1regular, Eigen::RowVector3d(0,0,0));
viewer.data().add_edges(PQC1regular, PQC2regular, Eigen::RowVector3d(0,0,0));
viewer.data().add_edges(PQC2regular, PQC3regular, Eigen::RowVector3d(0,0,0));
viewer.data().add_edges(PQC3regular, PQC0regular, Eigen::RowVector3d(0,0,0));
}
return false;
}
int main(int argc, char *argv[])
{
using namespace Eigen;
using namespace std;
// Load a quad mesh
igl::readOFF(TUTORIAL_SHARED_PATH "/halftunnel.off", VQC, FQC);
// Convert it in a triangle mesh
FQCtri.resize(2*FQC.rows(), 3);
FQCtri << FQC.col(0),FQC.col(1),FQC.col(2),
FQC.col(2),FQC.col(3),FQC.col(0);
igl::slice( VQC, FQC.col(0).eval(), 1, PQC0);
igl::slice( VQC, FQC.col(1).eval(), 1, PQC1);
igl::slice( VQC, FQC.col(2).eval(), 1, PQC2);
igl::slice( VQC, FQC.col(3).eval(), 1, PQC3);
// Create a planar version with ShapeUp
//igl::planarize_quad_mesh(VQC, FQC, 100, 0.005, VQCregular);
E.resize(FQC.size(),2);
E.col(0)<<FQC.col(0),FQC.col(1),FQC.col(2),FQC.col(3);
E.col(1)<<FQC.col(1),FQC.col(2),FQC.col(3),FQC.col(0);
VectorXi b(1); b(0)=0; //setting the first vertex to be the same.
VectorXd wShape=VectorXd::Constant(FQC.rows(),1.0);
VectorXd wSmooth=VectorXd::Constant(E.rows(),1.0);
MatrixXd bc(1,3); bc<<VQC.row(0);
VectorXi array_of_fours=VectorXi::Constant(FQC.rows(),4);
igl::shapeup_projection_function localFunction(igl::shapeup_regular_face_projection);
su_data.maxIterations=200;
shapeup_precomputation(VQC, array_of_fours,FQC,E,b,wShape, wSmooth,su_data);
shapeup_solve(bc,localFunction, VQC,su_data, false,VQCregular);
// Convert the planarized mesh to triangles
igl::slice( VQCregular, FQC.col(0).eval(), 1, PQC0regular);
igl::slice( VQCregular, FQC.col(1).eval(), 1, PQC1regular);
igl::slice( VQCregular, FQC.col(2).eval(), 1, PQC2regular);
igl::slice( VQCregular, FQC.col(3).eval(), 1, PQC3regular);
// Launch the viewer
igl::opengl::glfw::Viewer viewer;
key_down(viewer,'1',0);
viewer.data().invert_normals = true;
viewer.data().show_lines = false;
viewer.callback_key_down = &key_down;
viewer.launch();
}
例713）半隧道网格（左）已优化为几乎完全规则（右）。色标介于 [0,0.05] 之间，测量每个面的角度与 的平均归一化偏差 90∘
7.14 行进四面体 Marching Tetrahedra
- 通常，3D数据被捕获为在空间 f(x):R^3→R 上定义的标量场。潜伏在这个场中，标量场的等值面通常是突出的几何对象。值处的等值 v 面由 中的所有 R ^ 3 点 x 组成， f(x)=v 使得 .几何处理中的核心问题是将等值面提取为三角形网格，以便进一步进行基于网格的处理或可视化。这称为等轮廓。
- “行进四面体” 46 是在 3D 单纯复合体（又名 tet 网格）上等轮廓三线性函数 f 的著名方法。此方法的核心思想是使用从查找表中选择的预定义拓扑（也称为连通性）勾勒出通过每个像元（如果有的话）的等值面轮廓，具体取决于像元每个顶点处的函数值。该方法迭代（“行进”）复合体中的所有单元（“四面体”），并将最终网格拼接在一起。
- 在 libigl 中，构造一个三角形网格， igl::marching_tets 近似于在 tet 网格 (V,F) 位置 (TV, TT) 顶点采样的输入标量场 S isovalue 的值的等值水平集：
igl::marching_tets(TV,TT,S, isovalue ,V,F);
- 714案例完整展示：
#include <igl/opengl/glfw/Viewer.h>
#include <igl/readOBJ.h>
#include <igl/read_triangle_mesh.h>
#include <igl/voxel_grid.h>
#include <igl/tetrahedralized_grid.h>
#include <igl/marching_tets.h>
#include <igl/signed_distance.h>
#include <igl/writeOBJ.h>
#include <igl/get_seconds.h>
#include <igl/marching_cubes.h>
#include <Eigen/Core>
int main(int argc, char * argv[])
{
const auto & tictoc = []()
{
static double t_start = igl::get_seconds();
double diff = igl::get_seconds()-t_start;
t_start += diff;
return diff;
};
// Load a surface mesh which is a cube
Eigen::MatrixXd V;
Eigen::MatrixXi F;
igl::read_triangle_mesh(
argc>1?argv[1]: TUTORIAL_SHARED_PATH "/armadillo.obj",V,F);
Eigen::RowVector3i side;
Eigen::MatrixXd TV;
tictoc();
igl::voxel_grid(V,0,100,1,TV,side);
printf("igl::voxel_grid %g secs\n",tictoc());
Eigen::MatrixXi TT5,TT6;
tictoc();
igl::tetrahedralized_grid(TV,side,igl::TETRAHEDRALIZED_GRID_TYPE_5,TT5);
printf("igl::tetrahedralized_grid %g secs\n",tictoc());
igl::tetrahedralized_grid(TV,side,igl::TETRAHEDRALIZED_GRID_TYPE_6_ROTATIONAL,TT6);
printf("igl::tetrahedralized_grid %g secs\n",tictoc());
tictoc();
Eigen::VectorXd S;
{
Eigen::VectorXi I;
Eigen::MatrixXd C,N;
igl::signed_distance(
TV,V,F,
igl::SIGNED_DISTANCE_TYPE_FAST_WINDING_NUMBER,
std::numeric_limits<double>::min(),
std::numeric_limits<double>::max(),
S,I,C,N);
}
printf("igl::signed_distance %g secs\n",tictoc());
std::vector<Eigen::MatrixXd> SV(3);
std::vector<Eigen::MatrixXi> SF(3);
igl::marching_tets(TV,TT5,S,0,SV[0],SF[0]);
printf("igl::marching_tets5 %g secs\n",tictoc());
igl::marching_tets(TV,TT6,S,0,SV[1],SF[1]);
printf("igl::marching_tets6 %g secs\n",tictoc());
tictoc();
igl::marching_cubes(S,TV,side(0),side(1),side(2),0,SV[2],SF[2]);
printf("igl::marching_cubes %g secs\n",tictoc());
//igl::writeOBJ("tmc.obj",SV,SF);
// Draw the mesh stored in (SV, SF)
igl::opengl::glfw::Viewer vr;
vr.data().set_mesh(V,F);
vr.data().is_visible = false;
vr.append_mesh();
int sel = 0;
const auto update = [&]()
{
vr.data().clear();
vr.data().set_mesh(SV[sel],SF[sel]);
};
vr.callback_key_pressed = [&](decltype(vr) &,unsigned int key, int mod)
{
switch(key)
{
case ',':
case '.':
sel = (sel+SV.size()+(key=='.'?1:-1))%SV.size();
update();
return true;
case ' ':
vr.data_list[0].is_visible = !vr.data_list[0].is_visible;
vr.data_list[1].is_visible = !vr.data_list[1].is_visible;
vr.selected_data_index = (vr.selected_data_index+1)%2;
return true;
}
return false;
};
update();
vr.launch();
}
7.15 隐式函数网格划分 Implicit Function Meshing
- 715案例完整展示：
#include <igl/copyleft/marching_cubes.h>
#include <igl/unique_simplices.h>
#include <igl/dual_contouring.h>
#include <igl/get_seconds.h>
#include <igl/grid.h>
#include <igl/marching_cubes.h>
#include <igl/per_corner_normals.h>
#include <igl/parallel_for.h>
#include <igl/per_corner_normals.h>
#include <igl/per_face_normals.h>
#include <igl/polygon_corners.h>
#include <igl/slice.h>
#include <igl/sparse_voxel_grid.h>
#include <igl/opengl/glfw/Viewer.h>
int main(int argc, char * argv[])
{
const auto & tictoc = []()
{
static double t_start = igl::get_seconds();
double diff = igl::get_seconds()-t_start;
t_start += diff;
return diff;
};
// Create an interesting shape with sharp features using SDF CSG with spheres.
const auto & sphere = [](
const Eigen::RowVector3d & c,
const double r,
const Eigen::RowVector3d & x)->double
{
return (x-c).norm() - r;
};
const std::function<double(const Eigen::RowVector3d & x)> f =
[&](const Eigen::RowVector3d & x)->double
{
return
std::min(
std::min(std::max(
sphere(Eigen::RowVector3d(-0.2,0,-0.2),0.5,x),
-sphere(Eigen::RowVector3d(+0.2,0,0.2),0.5,x)),
sphere(Eigen::RowVector3d(-0.15,0,-0.15),0.3,x)
),
std::max(
std::max(
sphere(Eigen::RowVector3d(-0.2,-0.5,-0.2),0.6,x),x(1)+0.45),-0.6-x(1))
);
};
Eigen::RowVector3d p0(-0.2,0.5,-0.2);
assert(abs(f(p0)) < 1e-10 && "p0 should be on zero level-set");
// Simple finite difference gradients
const auto & fd = [](
const std::function<double(const Eigen::RowVector3d&)> & f,
const Eigen::RowVector3d & x)
{
const double eps = 1e-10;
Eigen::RowVector3d g;
for(int c = 0;c<3;c++)
{
const Eigen::RowVector3d xp = x+eps*Eigen::RowVector3d(c==0,c==1,c==2);
const double fp = f(xp);
const Eigen::RowVector3d xn = x-eps*Eigen::RowVector3d(c==0,c==1,c==2);
const double fn = f(xn);
g(c) = (fp-fn)/(2*eps);
}
return g;
};
const auto & f_grad = [&fd,&f](const Eigen::RowVector3d & x)
{
return fd(f,x).normalized();
};
Eigen::MatrixXd V;
Eigen::MatrixXi Q,F;
Eigen::MatrixXd mcV,mcN;
Eigen::MatrixXi mcF;
// Grid parameters
const Eigen::RowVector3d min_corner(-2,-2,-2);
const Eigen::RowVector3d max_corner(+2,+2,+2);
const int s = 256;
int nx = s+1;
int ny = s+1;
int nz = s+1;
const Eigen::RowVector3d step =
(max_corner-min_corner).array()/(Eigen::RowVector3d(nx,ny,nz).array()-1);
// Sparse grid below assumes regular grid
assert((step(0) == step(1))&&(step(0) == step(2)));
// Dual contouring parameters
bool constrained = false;
bool triangles = false;
bool root_finding = true;
for(int pass = 0;pass<2;pass++)
{
const bool sparse = pass == 1;
printf("Using %s grid..\n",sparse?"sparse":"dense");
if(sparse)
{
// igl::sparse_voxel_grid assumes (0,0,0) lies on the grid. But dense igl::grid
// below won't necessarily do that depending on nx,ny,nz.
tictoc();
Eigen::MatrixXd GV;
Eigen::VectorXd Gf;
Eigen::Matrix<int,Eigen::Dynamic,8> GI;
igl::sparse_voxel_grid(p0,f,step(0),16.*pow(step(0),-2.),Gf,GV,GI);
const auto t_Gf = tictoc();
printf(" %5f secs to populate sparse grid of %ld cells\n",t_Gf+tictoc(),GI.rows());
// Dual contouring requires list of sparse edges (not cells)
// extract _all_ edges from sparse_voxel_grid (conservative)
Eigen::Matrix<int,Eigen::Dynamic,2> GI2;
{
Eigen::Matrix<int,Eigen::Dynamic,2> all_GI2(GI.rows()*12,2);
all_GI2 <<
// front
GI.col(0),GI.col(1),
GI.col(1),GI.col(2),
GI.col(2),GI.col(3),
GI.col(3),GI.col(0),
// back
GI.col(4+0),GI.col(4+1),
GI.col(4+1),GI.col(4+2),
GI.col(4+2),GI.col(4+3),
GI.col(4+3),GI.col(4+0),
// sides
GI.col(0),GI.col(4+0),
GI.col(1),GI.col(4+1),
GI.col(2),GI.col(4+2),
GI.col(3),GI.col(4+3);
Eigen::VectorXi _1,_2;
igl::unique_simplices(all_GI2,GI2,_1,_2);
}
tictoc();
Eigen::RowVector3d step =
(max_corner-min_corner).array()/(Eigen::RowVector3d(nx,ny,nz).array()-1);
igl::dual_contouring(
f,f_grad,step,Gf,GV,GI2,constrained,triangles,root_finding,V,Q);
printf(" %5f secs dual contouring\n",t_Gf+tictoc());
// Could use igl::marching_cubes once
// https://github.com/libigl/libigl/pull/1687 is merged
tictoc();
igl::copyleft::marching_cubes(Gf,GV,GI,mcV, mcF);
printf(" %5f secs marching cubes\n",t_Gf+tictoc());
}else
{
tictoc();
igl::dual_contouring(
f,f_grad,min_corner,max_corner,nx,ny,nz,constrained,triangles,root_finding,V,Q);
printf(" %5f secs dual contouring\n",tictoc());
// build and sample grid
tictoc();
Eigen::MatrixXd GV;
igl::grid(Eigen::RowVector3i(nx,ny,nz),GV);
Eigen::VectorXd Gf(GV.rows());
igl::parallel_for(GV.rows(),[&](const int i)
{
GV.row(i).array() *= (max_corner-min_corner).array();
GV.row(i) += min_corner;
Gf(i) = f(GV.row(i));
},1000ul);
const auto t_grid = tictoc();
igl::marching_cubes(Gf,GV,nx,ny,nz,0,mcV,mcF);
const auto t_mc = tictoc();
printf(" %5f secs (%5f + %5f) marching cubes\n",t_grid+t_mc,t_grid,t_mc);
}
}
// Crisp (as possible) rendering of resulting MC triangle mesh
igl::per_corner_normals(mcV,mcF,20,mcN);
// Crisp rendering of resulting DC quad mesh with edges
Eigen::MatrixXi E;
Eigen::MatrixXd VV,N,NN;
Eigen::VectorXi J;
Eigen::MatrixXi FF;
if(triangles)
{
VV = V;
FF = Q;
E.resize(Q.rows()*3,2);
E<<
Q.col(0), Q.col(1),
Q.col(1), Q.col(2),
Q.col(2), Q.col(0);
}else
{
Eigen::VectorXi I,C;
igl::polygon_corners(Q,I,C);
E.resize(Q.rows()*4,2);
E<<
Q.col(0), Q.col(1),
Q.col(1), Q.col(2),
Q.col(2), Q.col(3),
Q.col(3), Q.col(0);
igl::per_face_normals(V,I,C,N,VV,FF,J);
igl::slice(N,J,1,NN);
igl::per_corner_normals(V,I,C,20,N,VV,FF,J,NN);
}
igl::opengl::glfw::Viewer vr;
bool show_edges = true;
bool use_dc = true;
const auto update = [&]()
{
const bool was_face_based = vr.data().face_based ;
vr.data().clear();
if(use_dc)
{
vr.data().set_mesh(VV,FF);
vr.data().show_lines = false;
vr.data().set_normals(NN);
if(show_edges)
{
vr.data().clear_edges();
vr.data().set_edges(V,E,Eigen::RowVector3d(0,0,0));
}
}else
{
vr.data().set_mesh(mcV,mcF);
vr.data().set_normals(mcN);
vr.data().show_lines = show_edges;
}
vr.data().face_based = was_face_based;
};
update();
vr.data().face_based = true;
vr.callback_key_pressed = [&](decltype(vr) &,unsigned int key, int mod)
{
switch(key)
{
case ' ': use_dc=!use_dc; update();return true;
case 'L': case 'l': show_edges=!show_edges; update();return true;
}
return false;
};
std::cout<<R"(
[space] Toggle between dual contouring and marching cubes
)";
vr.launch();
}
7.16 快速测地线距离近似的热法 Heat Method For Fast Geodesic Distance Approximation
-
在上面的精确离散测地线距离示例中，将精确计算测地线距离。这是一项昂贵的操作： O(n²log(n)) 对于带有边缘的 n 网格。2013年，Crane等人 47 提出了一种方法，通过求解表面上的热方程，过滤结果，然后通过求解泊松方程重建平滑解，更快地计算近似测地线距离。该方法从瓦拉丹的观察开始，两点 x y 之间的测地线距离 d(x,y) 等于一段时间 y 后扩散 x 的热量对数的平方根 t ：
-
热内核在哪里 kt,x 。我们可以将这种热扩散问题视为将热针放在上面 x ，然后在几秒钟后 t 测量点 y 的温度。
-
在三角形网格上，我们知道如何求解任何有限时间 t 的热方程：
-
如果我们有足够的数值精度和精度，我们可以简单地评估 √(−4tlogu) 一个小的时间参数 t 。Crane等人观察到的问题是，我们对数值 u 的数值精度远远不够。然而，渐变的方向 ∇u 出奇地准确。因此，他们的想法是获取 的 u 梯度，归一化这些向量以获得梯度方向（单位向量）。然后求解泊松方程，将这些方向积分到实际距离值中。
-
此方法涉及反转 n×n 稀疏矩阵（运算 O(n1.⋯) ），但如果使用 Cholesky 分解，则因式分解是预计算，即使更改了测地线距离的来源，也可以重用。对于新源，只需执行反向替换。
-
在 libigl 中，您可以使用此方法通过两个步骤计算网格 （ V ， F ） 从源顶点索引列表 gamma 到矢量 D 的近似测地线距离：
igl::HeatGeodesicsData<double> data;
igl::heat_geodesics_precompute(V,F,data);
...
igl::heat_geodesics_solve(data,gamma,D);
- 715案例完整展示：
#include <igl/read_triangle_mesh.h>
#include <igl/triangulated_grid.h>
#include <igl/heat_geodesics.h>
#include <igl/unproject_onto_mesh.h>
#include <igl/avg_edge_length.h>
#include <igl/isolines_map.h>
#include <igl/opengl/glfw/Viewer.h>
#include <igl/opengl/create_shader_program.h>
#include <igl/opengl/destroy_shader_program.h>
#include <iostream>
void set_colormap(igl::opengl::glfw::Viewer & viewer)
{
const int num_intervals = 30;
Eigen::MatrixXd CM(num_intervals,3);
// Colormap texture
for(int i = 0;i<num_intervals;i++)
{
double t = double(num_intervals - i - 1)/double(num_intervals-1);
CM(i,0) = std::max(std::min(2.0*t-0.0,1.0),0.0);
CM(i,1) = std::max(std::min(2.0*t-1.0,1.0),0.0);
CM(i,2) = std::max(std::min(6.0*t-5.0,1.0),0.0);
}
igl::isolines_map(Eigen::MatrixXd(CM),CM);
viewer.data().set_colormap(CM);
}
int main(int argc, char *argv[])
{
// Create the peak height field
Eigen::MatrixXi F;
Eigen::MatrixXd V;
igl::read_triangle_mesh( argc>1?argv[1]: TUTORIAL_SHARED_PATH "/beetle.off",V,F);
// Precomputation
igl::HeatGeodesicsData<double> data;
double t = std::pow(igl::avg_edge_length(V,F),2);
const auto precompute = [&]()
{
if(!igl::heat_geodesics_precompute(V,F,t,data))
{
std::cerr<<"Error: heat_geodesics_precompute failed."<<std::endl;
exit(EXIT_FAILURE);
};
};
precompute();
igl::opengl::glfw::Viewer viewer;
bool down_on_mesh = false;
const auto update = [&]()->bool
{
int fid;
Eigen::Vector3f bc;
// Cast a ray in the view direction starting from the mouse position
double x = viewer.current_mouse_x;
double y = viewer.core().viewport(3) - viewer.current_mouse_y;
if(igl::unproject_onto_mesh(Eigen::Vector2f(x,y), viewer.core().view,
viewer.core().proj, viewer.core().viewport, V, F, fid, bc))
{
Eigen::VectorXd D;
// if big mesh, just use closest vertex. Otherwise, blend distances to
// vertices of face using barycentric coordinates.
if(F.rows()>100000)
{
// 3d position of hit
const Eigen::RowVector3d m3 =
V.row(F(fid,0))*bc(0) + V.row(F(fid,1))*bc(1) + V.row(F(fid,2))*bc(2);
int cid = 0;
Eigen::Vector3d(
(V.row(F(fid,0))-m3).squaredNorm(),
(V.row(F(fid,1))-m3).squaredNorm(),
(V.row(F(fid,2))-m3).squaredNorm()).minCoeff(&cid);
const int vid = F(fid,cid);
igl::heat_geodesics_solve(data,(Eigen::VectorXi(1,1)<<vid).finished(),D);
}else
{
D = Eigen::VectorXd::Zero(V.rows());
for(int cid = 0;cid<3;cid++)
{
const int vid = F(fid,cid);
Eigen::VectorXd Dc;
igl::heat_geodesics_solve(data,(Eigen::VectorXi(1,1)<<vid).finished(),Dc);
D += Dc*bc(cid);
}
}
viewer.data().set_data(D);
return true;
}
return false;
};
viewer.callback_mouse_down =
[&](igl::opengl::glfw::Viewer& viewer, int, int)->bool
{
if(update())
{
down_on_mesh = true;
return true;
}
return false;
};
viewer.callback_mouse_move =
[&](igl::opengl::glfw::Viewer& viewer, int, int)->bool
{
if(down_on_mesh)
{
update();
return true;
}
return false;
};
viewer.callback_mouse_up =
[&down_on_mesh](igl::opengl::glfw::Viewer& viewer, int, int)->bool
{
down_on_mesh = false;
return false;
};
std::cout<<R"(Usage:
[click] Click on shape to pick new geodesic distance source
,/. Decrease/increase t by factor of 10.0
D,d Toggle using intrinsic Delaunay discrete differential operators
)";
viewer.callback_key_pressed =
[&](igl::opengl::glfw::Viewer& /*viewer*/, unsigned int key, int mod)->bool
{
switch(key)
{
default:
return false;
case 'D':
case 'd':
data.use_intrinsic_delaunay = !data.use_intrinsic_delaunay;
std::cout<<(data.use_intrinsic_delaunay?"":"not ")<<
"using intrinsic delaunay..."<<std::endl;
precompute();
update();
break;
case '.':
case ',':
t *= (key=='.'?10.0:0.1);
precompute();
update();
std::cout<<"t: "<<t<<std::endl;
break;
}
return true;
};
// Show mesh
viewer.data().set_mesh(V, F);
viewer.data().set_data(Eigen::VectorXd::Zero(V.rows()));
set_colormap(viewer);
viewer.data().show_lines = false;
viewer.launch();
}
加载网格并计算用户单击位置的近似测地线距离
7.16.1 内在德劳奈三角测量 Intrinsic Delaunay Triangulation
- 测地线距离的原始热法在规则的无偏网格上效果很好：即，有限元余切矩阵和质量矩阵表现良好。但是，对于质量较差的网格，此方法可能会显示任意差的结果。增加时间参数 t 可以减少这种不稳定性，但同时可以平滑生成的近似距离。
- 相反，改进的一种途径是使用所谓的内在 Delaunay 三角测量离散拉普拉斯算子 48 。
- 由于余切拉普拉斯仅取决于三角形网格的边长，因此此新算子将通过固有地翻转边并记录边长度的变化来构造。翻转边，直到每条边都是局部 Delaunay （即，其相应的余切权重为正）。
- 您可以使用以下命令计算网格 （ 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);
7.17 汤和云的快速缠绕数 Fast Winding Number For Soups And Clouds
2018年，Barill等人 50 演示了如何显著加快上述广义绕组数的计算速度。三角形网格的广义缠绕数的原始定义也扩展到（定向）点云。
7.17.1 汤 Soups
- 对于三角形汤， 40 理想的精确分而治之的方法在三角形的数量上以对数缩放 O(logn) 。但是，此方法具有计算变为线性 O(n) 的故障模式。使用类似于n体系统模拟或静电问题中使用的树方法可以显着改善这一点。结果是近似的，但更紧密地遵循趋势 O(logn) ，并且常数因子要小得多。
- 计算汤的快速缠绕数有两个步骤：构建树形数据结构，然后在查询点进行评估。在 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);
- 717案例完整展示：
#include <igl/fast_winding_number.h>
#include <igl/read_triangle_mesh.h>
#include <igl/slice_mask.h>
#include <Eigen/Geometry>
#include <igl/octree.h>
#include <igl/barycenter.h>
#include <igl/knn.h>
#include <igl/random_points_on_mesh.h>
#include <igl/bounding_box_diagonal.h>
#include <igl/per_face_normals.h>
#include <igl/copyleft/cgal/point_areas.h>
#include <igl/opengl/glfw/Viewer.h>
#include <igl/get_seconds.h>
#include <iostream>
#include <cstdlib>
int main(int argc, char *argv[])
{
const auto time = [](std::function<void(void)> func)->double
{
const double t_before = igl::get_seconds();
func();
const double t_after = igl::get_seconds();
return t_after-t_before;
};
Eigen::MatrixXd V;
Eigen::MatrixXi F;
igl::read_triangle_mesh(argc>1?argv[1]:TUTORIAL_SHARED_PATH "/bunny.off",V,F);
// Sample mesh for point cloud
Eigen::MatrixXd P,N;
{
Eigen::VectorXi I;
Eigen::SparseMatrix<double> B;
igl::random_points_on_mesh(10000,V,F,B,I);
P = B*V;
Eigen::MatrixXd FN;
igl::per_face_normals(V,F,FN);
N.resize(P.rows(),3);
for(int p = 0;p<I.rows();p++)
{
N.row(p) = FN.row(I(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);
}
if(argc<=1)
{
// corrupt mesh
Eigen::MatrixXd BC;
igl::barycenter(V,F,BC);
Eigen::MatrixXd OV = V;
V.resize(F.rows()*3,3);
for(int f = 0;f<F.rows();f++)
{
for(int c = 0;c<3;c++)
{
int v = f+c*F.rows();
// random rotation about barycenter
Eigen::AngleAxisd R(
0.5*static_cast <double> (rand()) / static_cast <double> (RAND_MAX),
Eigen::Vector3d::Random(3,1));
V.row(v) = (OV.row(F(f,c))-BC.row(f))*R.matrix()+BC.row(f);
F(f,c) = v;
}
}
}
// Generate a list of random query points in the bounding box
Eigen::MatrixXd Q = Eigen::MatrixXd::Random(1000000,3);
const Eigen::RowVector3d Vmin = V.colwise().minCoeff();
const Eigen::RowVector3d Vmax = V.colwise().maxCoeff();
const Eigen::RowVector3d Vdiag = Vmax-Vmin;
for(int q = 0;q<Q.rows();q++)
{
Q.row(q) = (Q.row(q).array()*0.5+0.5)*Vdiag.array() + Vmin.array();
}
// Positions of points inside of point cloud P
Eigen::MatrixXd QiP;
{
Eigen::MatrixXd O_CM;
Eigen::VectorXd O_R;
Eigen::MatrixXd O_EC;
printf(" point cloud precomputation (% 8ld points): %g secs\n",
P.rows(),
time([&](){igl::fast_winding_number(P,N,A,O_PI,O_CH,2,O_CM,O_R,O_EC);}));
Eigen::VectorXd WiP;
printf(" point cloud evaluation (% 8ld queries): %g secs\n",
Q.rows(),
time([&](){igl::fast_winding_number(P,N,A,O_PI,O_CH,O_CM,O_R,O_EC,Q,2,WiP);}));
igl::slice_mask(Q,(WiP.array()>0.5).eval(),1,QiP);
}
// Positions of points inside of triangle soup (V,F)
Eigen::MatrixXd QiV;
{
igl::FastWindingNumberBVH fwn_bvh;
printf("triangle soup precomputation (% 8ld triangles): %g secs\n",
F.rows(),
time([&](){igl::fast_winding_number(V.cast<float>().eval(),F,2,fwn_bvh);}));
Eigen::VectorXf WiV;
printf(" triangle soup evaluation (% 8ld queries): %g secs\n",
Q.rows(),
time([&](){igl::fast_winding_number(fwn_bvh,2,Q.cast<float>().eval(),WiV);}));
igl::slice_mask(Q,WiV.array()>0.5,1,QiV);
}
// Visualization
igl::opengl::glfw::Viewer viewer;
// For dislpaying normals as little line segments
Eigen::MatrixXd PN(2*P.rows(),3);
Eigen::MatrixXi E(P.rows(),2);
const double bbd = igl::bounding_box_diagonal(V);
for(int p = 0;p<P.rows();p++)
{
E(p,0) = 2*p;
E(p,1) = 2*p+1;
PN.row(E(p,0)) = P.row(p);
PN.row(E(p,1)) = P.row(p)+bbd*0.01*N.row(p);
}
bool show_P = false;
int show_Q = 0;
int query_data = 0;
viewer.data_list[query_data].set_mesh(V,F);
viewer.data_list[query_data].clear();
viewer.data_list[query_data].point_size = 2;
viewer.append_mesh();
int object_data = 1;
viewer.data_list[object_data].set_mesh(V,F);
viewer.data_list[object_data].point_size = 5;
const auto update = [&]()
{
viewer.data_list[query_data].clear();
switch(show_Q)
{
case 1:
// show all Q
viewer.data_list[query_data].set_points(Q,Eigen::RowVector3d(0.996078,0.760784,0.760784));
break;
case 2:
// show all Q inside
if(show_P)
{
viewer.data_list[query_data].set_points(QiP,Eigen::RowVector3d(0.564706,0.847059,0.768627));
}else
{
viewer.data_list[query_data].set_points(QiV,Eigen::RowVector3d(0.564706,0.847059,0.768627));
}
break;
}
viewer.data_list[object_data].clear();
if(show_P)
{
viewer.data_list[object_data].set_points(P,Eigen::RowVector3d(1,1,1));
viewer.data_list[object_data].set_edges(PN,E,Eigen::RowVector3d(0.8,0.8,0.8));
}else
{
viewer.data_list[object_data].set_mesh(V,F);
}
};
viewer.callback_key_pressed =
[&](igl::opengl::glfw::Viewer &, unsigned int key, int mod)
{
switch(key)
{
default:
return false;
case '1':
show_P = !show_P;
break;
case '2':
show_Q = (show_Q+1) % 3;
break;
}
update();
return true;
};
std::cout<<R"(
FastWindingNumber
1 Toggle point cloud and triangle soup
2 Toggle hiding query points, showing query points, showing inside queries
)";
update();
viewer.launch();
}
- 键盘1
- 键盘2
加载一个网格，对 1000000 个随机查询进行采样，然后丢弃给定模型之外的所有查询。（图片/兔子汤.jpg）
7.17.2 云 Clouds
- 对于点云，过程类似，但可能需要更多的预计算。定义点云的缠绕数，只要每个点带有一个朝外的法线和一个面积值。可以使用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);
7.18 迭代最近点 Iterative Closest Point¶
- 718案例完整展示：
// This file is part of libigl, a simple c++ geometry processing library.
//
// Copyright (C) 2019 Alec Jacobson <alecjacobson@gmail.com>
//
// This Source Code Form is subject to the terms of the Mozilla Public License
// v. 2.0. If a copy of the MPL was not distributed with this file, You can
// obtain one at http://mozilla.org/MPL/2.0/.
#include <Eigen/Core>
#include <Eigen/Geometry>
#include <igl/read_triangle_mesh.h>
#include <igl/opengl/glfw/Viewer.h>
#include <igl/iterative_closest_point.h>
#include <igl/random_dir.h>
#include <igl/PI.h>
#include <igl/AABB.h>
#include <igl/per_face_normals.h>
#include <igl/avg_edge_length.h>
#include <igl/per_vertex_normals.h>
#include <igl/barycenter.h>
#include <igl/slice.h>
#include <igl/slice_mask.h>
#include <iostream>
int main(int argc, char * argv[])
{
Eigen::MatrixXd OVX,VX,VY;
Eigen::MatrixXi FX,FY;
igl::read_triangle_mesh( argc>1?argv[1]: TUTORIAL_SHARED_PATH "/decimated-max.obj",VY,FY);
const double bbd = (VY.colwise().maxCoeff()-VY.colwise().minCoeff()).norm();
FX = FY;
{
// sprinkle a noise so that we can see z-fighting when the match is perfect.
const double h = igl::avg_edge_length(VY,FY);
OVX = VY + 1e-2*h*Eigen::MatrixXd::Random(VY.rows(),VY.cols());
}
VX = OVX;
igl::AABB<Eigen::MatrixXd,3> Ytree;
Ytree.init(VY,FY);
Eigen::MatrixXd NY;
igl::per_face_normals(VY,FY,NY);
igl::opengl::glfw::Viewer v;
std::cout<<R"(
[space] conduct a single iterative closest point iteration
R,r reset to a random orientation and offset
)";
const auto apply_random_rotation = [&]()
{
const Eigen::Matrix3d R = Eigen::AngleAxisd(
2.*igl::PI*(double)rand()/RAND_MAX*0.3, igl::random_dir()).matrix();
const Eigen::RowVector3d cen =
0.5*(VY.colwise().maxCoeff()+VY.colwise().minCoeff());
VX = ((OVX*R).rowwise()+(cen-cen*R)).eval();
};
const auto single_iteration = [&]()
{
// Perform single iteration of ICP method
Eigen::Matrix3d R;
Eigen::RowVector3d t;
igl::iterative_closest_point(VX,FX,VY,FY,Ytree,NY,1000,1,R,t);
VX = ((VX*R).rowwise()+t).eval();
v.data().set_mesh(VX,FX);
v.data().compute_normals();
};
v.callback_pre_draw = [&](igl::opengl::glfw::Viewer &)->bool
{
if(v.core().is_animating)
{
single_iteration();
}
return false;
};
v.callback_key_pressed =
[&](igl::opengl::glfw::Viewer &,unsigned char key,int)->bool
{
switch(key)
{
case ' ':
{
v.core().is_animating = false;
single_iteration();
return true;
}
case 'R':
case 'r':
// Random rigid transformation
apply_random_rotation();
v.data().set_mesh(VX,FX);
v.data().compute_normals();
return true;
break;
}
return false;
};
v.data().set_mesh(VY,FY);
v.data().set_colors(Eigen::RowVector3d(1,1,1));
v.data().show_lines = false;
v.append_mesh();
v.data().set_mesh(VX,FX);
v.data().show_lines = false;
v.launch();
}
7.19 分解图 Exploded View
- 719案例完整展示：
// This file is part of libigl, a simple c++ geometry processing library.
//
// Copyright (C) 2020 Alec Jacobson <alecjacobson@gmail.com>
//
// This Source Code Form is subject to the terms of the Mozilla Public License
// v. 2.0. If a copy of the MPL was not distributed with this file, You can
// obtain one at http://mozilla.org/MPL/2.0/.
#include <igl/opengl/glfw/Viewer.h>
#include <igl/readMESH.h>
#include <igl/barycenter.h>
#include <igl/volume.h>
#include <igl/exploded_view.h>
#include <igl/slice.h>
#include <igl/colormap.h>
int main(int argc, char *argv[])
{
Eigen::MatrixXd V;
Eigen::MatrixXi T,F;
igl::readMESH(argc>1?argv[1]: TUTORIAL_SHARED_PATH "/octopus-low.mesh",V,T,F);
// Some per-tet data
Eigen::VectorXd D;
{
Eigen::MatrixXd BC;
igl::barycenter(V,T,BC);
Eigen::VectorXd vol;
igl::volume(V,T,vol);
const Eigen::RowVectorXd c = vol.transpose()*BC/vol.array().sum();
D = (BC.rowwise()-c).rowwise().norm();
}
// Plot the mesh
igl::opengl::glfw::Viewer viewer;
double t = 1;
double s = 1;
const auto update = [&]()
{
Eigen::MatrixXd EV;
Eigen::MatrixXi EF;
Eigen::VectorXi I,J;
igl::exploded_view(V,T,s,t,EV,EF,I,J);
Eigen::VectorXd DJ;
igl::slice(D,J,1,DJ);
static bool first = true;
if(first)
{
viewer.data().clear();
viewer.data().set_mesh(EV,EF);
first = false;
}else
{
viewer.data().set_vertices(EV);
}
viewer.data().set_face_based(true);
Eigen::MatrixXd C;
igl::colormap(igl::COLOR_MAP_TYPE_VIRIDIS,DJ,true,C);
viewer.data().set_colors(C);
};
int mod = 0;
float prev_y;
viewer.callback_mouse_move = [&](igl::opengl::glfw::Viewer &, int x, int y)->bool
{
if((mod & IGL_MOD_SHIFT)||(mod & IGL_MOD_ALT))
{
if(mod & IGL_MOD_SHIFT)
{
t = std::min(std::max(t+0.001*(y-prev_y),1.),2.);
}
if(mod & IGL_MOD_ALT)
{
s = std::min(std::max(s+0.001*(y-prev_y),0.),1.);
}
prev_y = y;
update();
return true;
}
return false;
};
// get modifier
viewer.callback_key_down =
[&](igl::opengl::glfw::Viewer &, unsigned char key, int _mod)->bool
{
prev_y = viewer.current_mouse_y;
mod = _mod;
return false;
};
viewer.callback_key_up =
[&](igl::opengl::glfw::Viewer &, unsigned char key, int _mod)->bool
{
mod = _mod;
return false;
};
update();
viewer.launch();
}
7.20 蓝色噪声表面采样 Blue Noise Surface Sampling
- 720案例完整展示：
// This file is part of libigl, a simple c++ geometry processing library.
//
// Copyright (C) 2020 Alec Jacobson <alecjacobson@gmail.com>
//
// This Source Code Form is subject to the terms of the Mozilla Public License
// v. 2.0. If a copy of the MPL was not distributed with this file, You can
// obtain one at http://mozilla.org/MPL/2.0/.
#include <igl/opengl/glfw/Viewer.h>
#include <igl/read_triangle_mesh.h>
#include <igl/blue_noise.h>
#include <igl/per_vertex_normals.h>
#include <igl/barycentric_interpolation.h>
#include <igl/blue_noise.h>
#include <igl/doublearea.h>
#include <igl/random_points_on_mesh.h>
#include <igl/PI.h>
int main(int argc, char *argv[])
{
Eigen::MatrixXd V;
Eigen::MatrixXi F;
igl::read_triangle_mesh(
argc>1?argv[1]: TUTORIAL_SHARED_PATH "/elephant.obj",V,F);
Eigen::MatrixXd N;
igl::per_vertex_normals(V,F,N);
const double bbd = (V.colwise().maxCoeff()- V.colwise().minCoeff()).norm();
Eigen::MatrixXd P_blue;
Eigen::MatrixXd N_blue;
Eigen::VectorXd A_blue;
Eigen::MatrixXd C_blue;
{
const int n_desired = argc>2?atoi(argv[2]):50000;
// Heuristic to determine radius from desired number
const double r = [&V,&F](const int n)
{
Eigen::VectorXd A;
igl::doublearea(V,F,A);
return sqrt(((A.sum()*0.5/(n*0.6162910373))/igl::PI));
}(n_desired);
printf("blue noise radius: %g\n",r);
Eigen::MatrixXd B;
Eigen::VectorXi I;
igl::blue_noise(V,F,r,B,I,P_blue);
igl::barycentric_interpolation(N,F,B,I,N_blue);
N_blue.rowwise().normalize();
}
Eigen::MatrixXd P_white;
Eigen::MatrixXd N_white;
Eigen::VectorXd A_white;
Eigen::MatrixXd C_white;
{
Eigen::MatrixXd B;
Eigen::VectorXi I;
igl::random_points_on_mesh(P_blue.rows(),V,F,B,I,P_white);
igl::barycentric_interpolation(N,F,B,I,N_white);
N_white.rowwise().normalize();
}
// Color
Eigen::RowVector3d C(1,1,1);
// Plot the mesh
igl::opengl::glfw::Viewer viewer;
viewer.data().set_mesh(V,F);
viewer.data().show_lines = false;
viewer.data().show_faces = false;
viewer.data().point_size = 4;
bool use_blue = true;
const auto update = [&]()
{
const Eigen::RowVector3d orange(1.0,0.7,0.2);
if(use_blue)
{
viewer.data().set_points(P_blue,Eigen::RowVector3d(0.3,0.4,1.0));
viewer.data().set_edges_from_vector_field(P_blue,bbd*0.01*N_blue,orange);
}else
{
viewer.data().set_points(P_white,Eigen::RowVector3d(1,1,1));
viewer.data().set_edges_from_vector_field(P_white,bbd*0.01*N_white,orange);
}
};
update();
viewer.callback_key_pressed =
[&](igl::opengl::glfw::Viewer &,unsigned char key,int)->bool
{
switch(key)
{
case ' ':
{
use_blue = !use_blue;
update();
return true;
}
}
return false;
};
std::cout<<R"(
[space] toggle between blue and white noise samplings
)";
viewer.launch();
}
7.21 矢量场平滑 Vector Field Smoothing
- 721案例完整展示：
#include <igl/read_triangle_mesh.h>
#include <igl/parula.h>
#include <igl/remove_unreferenced.h>
#include <igl/opengl/glfw/Viewer.h>
#include <igl/per_face_normals.h>
#include <igl/orient_halfedges.h>
#include <igl/cr_vector_laplacian.h>
#include <igl/cr_vector_mass.h>
#include <igl/edge_midpoints.h>
#include <igl/edge_vectors.h>
#include <igl/average_from_edges_onto_vertices.h>
#include <igl/PI.h>
#include <Eigen/Core>
#include <Eigen/SparseCholesky>
#include <Eigen/Geometry>
#include <iostream>
#include <set>
#include <limits>
#include <stdlib.h>
int main(int argc, char * argv[])
{
typedef Eigen::SparseMatrix<double> SparseMat;
//Constants used for smoothing
const double howMuchToSmoothBy = 1e-1;
const int howManySmoothingInterations = 50;
//Read our mesh
Eigen::MatrixXd V;
Eigen::MatrixXi F;
if(!igl::read_triangle_mesh
(argc>1?argv[1]: TUTORIAL_SHARED_PATH "/elephant.obj",V,F)) {
std::cout << "Failed to load mesh." << std::endl;
}
//Orient edges for plotting
Eigen::MatrixXi E, oE;
igl::orient_halfedges(F, E, oE);
//Compute edge midpoints & edge vectors
Eigen::MatrixXd edgeMps, parVec, perpVec;
igl::edge_midpoints(V, F, E, oE, edgeMps);
igl::edge_vectors(V, F, E, oE, parVec, perpVec);
//Constructing a function to add noise to
const auto zraw_function = [] (const Eigen::Vector3d& x) {
return Eigen::Vector3d(0.2*x(1) + cos(2*x(1)+0.2),
0.5*x(0) + 0.15,
0.3*cos(0.2+igl::PI*x(2)));
};
Eigen::VectorXd zraw(2*edgeMps.rows());
for(int i=0; i<edgeMps.rows(); ++i) {
const Eigen::Vector3d f = zraw_function(edgeMps.row(i));
zraw(i) = f.dot(parVec.row(i));
zraw(i+edgeMps.rows()) = f.dot(perpVec.row(i));
}
//Add noise
srand(71);
const double l = 15;
Eigen::VectorXd znoisy = zraw + l*Eigen::VectorXd::Random(zraw.size());
//Denoise function using the vector Dirichlet energy
Eigen::VectorXd zsmoothed = znoisy;
for(int i=0; i<howManySmoothingInterations; ++i) {
//Compute Laplacian and mass matrix
SparseMat L, M;
igl::cr_vector_mass(V, F, E, oE, M);
igl::cr_vector_laplacian(V, F, E, oE, L);
//Implicit step
Eigen::SimplicialLDLT<SparseMat> rhsSolver(M + howMuchToSmoothBy*L);
zsmoothed = rhsSolver.solve(M*zsmoothed);
}
//Convert vector fields for plotting
const auto cr_result_to_vecs_and_colors = [&]
(const Eigen::VectorXd& z, Eigen::MatrixXd& vecs, Eigen::MatrixXd& colors) {
vecs.resize(edgeMps.rows(), 3);
for(int i=0; i<edgeMps.rows(); ++i) {
vecs.row(i) = z(i)*parVec.row(i)
+ z(i+edgeMps.rows())*perpVec.row(i);
}
igl::average_from_edges_onto_vertices
(F, E, oE, vecs.rowwise().norm(), colors);
};
Eigen::MatrixXd noisyvecs, noisycolors, smoothedvecs, smoothedcolors,
rawvecs, rawcolors;
cr_result_to_vecs_and_colors(znoisy, noisyvecs, noisycolors);
cr_result_to_vecs_and_colors(zsmoothed, smoothedvecs, smoothedcolors);
cr_result_to_vecs_and_colors(zraw, rawvecs, rawcolors);
//Viewer that shows noisy and denoised functions
igl::opengl::glfw::Viewer viewer;
viewer.data().set_mesh(V,F);
viewer.data().show_lines = false;
viewer.callback_key_down =
[&](igl::opengl::glfw::Viewer & viewer, unsigned char key, int mod)->bool
{
const Eigen::MatrixXd *vecs, *colors;
switch(key) {
case '1':
vecs = &rawvecs;
colors = &rawcolors;
break;
case '2':
vecs = &noisyvecs;
colors = &noisycolors;
break;
case '3':
vecs = &smoothedvecs;
colors = &smoothedcolors;
break;
default:
return false;
}
viewer.data().set_data(*colors);
viewer.data().clear_edges();
const double s = 0.08; //How much to scale vectors during plotting
viewer.data().add_edges(edgeMps, edgeMps + s*(*vecs),
Eigen::RowVector3d(0.1, 0.1, 0.1));
return true;
};
std::cout << R"(Usage:
1 Show raw function
2 Show noisy function
3 Show smoothed function
)";
Eigen::MatrixXd CM;
igl::parula(Eigen::VectorXd::LinSpaced
(500,znoisy.minCoeff(), znoisy.maxCoeff()).eval(), true, CM);
viewer.data().set_colormap(CM);
viewer.callback_key_down(viewer, '1', 0);
viewer.launch();
return 0;
}
7.22 向量并行传输 Vector Parallel Transport
- 722案例完整展示：
#include <igl/read_triangle_mesh.h>
#include <igl/parula.h>
#include <igl/remove_unreferenced.h>
#include <igl/opengl/glfw/Viewer.h>
#include <igl/per_face_normals.h>
#include <igl/orient_halfedges.h>
#include <igl/cr_vector_laplacian.h>
#include <igl/cr_vector_mass.h>
#include <igl/crouzeix_raviart_cotmatrix.h>
#include <igl/crouzeix_raviart_massmatrix.h>
#include <igl/edge_midpoints.h>
#include <igl/edge_vectors.h>
#include <igl/average_from_edges_onto_vertices.h>
#include <igl/min_quad_with_fixed.h>
#include <igl/heat_geodesics.h>
#include <Eigen/Core>
#include <Eigen/SparseCholesky>
#include <Eigen/Geometry>
#include <iostream>
#include <set>
#include <limits>
#include <stdlib.h>
int main(int argc, char * argv[])
{
typedef Eigen::SparseMatrix<double> SparseMat;
typedef Eigen::Matrix<double, 1, 1> Vector1d;
typedef Eigen::Matrix<int, 1, 1> Vector1i;
//Constants used for smoothing
const double howMuchToSmoothBy = 1e-1;
const int howManySmoothingInterations = 50;
//Read our mesh
Eigen::MatrixXd V;
Eigen::MatrixXi F;
if(!igl::read_triangle_mesh
(argc>1?argv[1]: TUTORIAL_SHARED_PATH "/cheburashka.off",V,F)) {
std::cout << "Failed to load mesh." << std::endl;
}
//Compute vector Laplacian and mass matrix
Eigen::MatrixXi E, oE;//Compute Laplacian and mass matrix
SparseMat vecL, vecM;
igl::cr_vector_mass(V, F, E, oE, vecM);
igl::cr_vector_laplacian(V, F, E, oE, vecL);
const int m = vecL.rows()/2; //The number of edges in the mesh
//Convert the E / oE matrix format to list of edges / EMAP format required
// by the functions constructing scalar Crouzeix-Raviart functions
Eigen::MatrixXi Elist(m,2), EMAP(3*F.rows(),1);
for(int i=0; i<F.rows(); ++i) {
for(int j=0; j<3; ++j) {
const int e = E(i,j);
EMAP(i+j*F.rows()) = e;
if(oE(i,j)>0) {
Elist.row(e) << F(i, (j+1)%3), F(i, (j+2)%3);
}
}
}
SparseMat scalarL, scalarM;
igl::crouzeix_raviart_massmatrix(V, F, Elist, EMAP, scalarM);
igl::crouzeix_raviart_cotmatrix(V, F, Elist, EMAP, scalarL);
//Compute edge midpoints & edge vectors
Eigen::MatrixXd edgeMps, parVec, perpVec;
igl::edge_midpoints(V, F, E, oE, edgeMps);
igl::edge_vectors(V, F, E, oE, parVec, perpVec);
//Perform the vector heat method
const int initialIndex = 14319;
const double initialPara=0.95, initialPerp=0.08;
const double t = 0.01;
SparseMat Aeq;
Eigen::VectorXd Beq;
Eigen::VectorXi known = Eigen::Vector2i(initialIndex, initialIndex+m);
Eigen::VectorXd knownVals = Eigen::Vector2d(initialPara, initialPerp);
Eigen::VectorXd Y0 = Eigen::VectorXd::Zero(2*m), Yt;
Y0(initialIndex) = initialPara; Y0(initialIndex+m) = initialPerp;
igl::min_quad_with_fixed
(SparseMat(vecM+t*vecL), Eigen::VectorXd(-vecM*Y0), known, knownVals,
Aeq, Beq, false, Yt);
Eigen::VectorXd u0 = Eigen::VectorXd::Zero(m), ut;
u0(initialIndex) = sqrt(initialPara*initialPara + initialPerp*initialPerp);
Eigen::VectorXi knownScal = Vector1i(initialIndex);
Eigen::VectorXd knownScalVals = Vector1d(u0(initialIndex));
igl::min_quad_with_fixed
(SparseMat(scalarM+t*scalarL), Eigen::VectorXd(-scalarM*u0), knownScal,
knownScalVals, Aeq, Beq, false, ut);
Eigen::VectorXd phi0 = Eigen::VectorXd::Zero(m), phit;
phi0(initialIndex) = 1;
Eigen::VectorXd knownScalValsPhi = Vector1d(1);
igl::min_quad_with_fixed
(SparseMat(scalarM+t*scalarL), Eigen::VectorXd(-scalarM*phi0), knownScal,
knownScalValsPhi, Aeq, Beq, false, phit);
Eigen::ArrayXd Xtfactor = ut.array() /
(phit.array() * (Yt.array().segment(0,m)*Yt.array().segment(0,m)
+ Yt.array().segment(m,m)*Yt.array().segment(m,m)).sqrt());
Eigen::VectorXd Xt(2*m);
Xt.segment(0,m) = Xtfactor * Yt.segment(0,m).array();
Xt.segment(m,m) = Xtfactor * Yt.segment(m,m).array();
//Compute scalar heat colors
igl::HeatGeodesicsData<double> hgData;
igl::heat_geodesics_precompute(V, F, hgData);
Eigen::VectorXd heatColor;
Eigen::VectorXi gamma = Elist.row(initialIndex);
igl::heat_geodesics_solve(hgData, gamma, heatColor);
//Convert vector field for plotting
Eigen::MatrixXd vecs(m, 3);
for(int i=0; i<edgeMps.rows(); ++i) {
vecs.row(i) = Xt(i)*parVec.row(i) + Xt(i+edgeMps.rows())*perpVec.row(i);
}
//Viewer that shows parallel transported vector
igl::opengl::glfw::Viewer viewer;
viewer.data().set_mesh(V,F);
viewer.data().show_lines = false;
viewer.data().set_data(heatColor.maxCoeff()-heatColor.array(), //invert colormap
igl::COLOR_MAP_TYPE_VIRIDIS);
const double s = 0.012; //How much to scale vectors during plotting
Eigen::MatrixXd vecColors(m, 3);
for(int i=0; i<m; ++i) {
vecColors.row(i) << 0.1, 0.1, 0.1;
}
vecColors.row(initialIndex) << 0.9, 0.1, 0.1;
viewer.data().add_edges(edgeMps, edgeMps + s*vecs, vecColors);
std::cout << R"(The red vector is parallel transported to every point on the surface.
The surface is shaded by geodesic distance from the red vector.
)"
<< std::endl;
viewer.launch();
return 0;
}
5454
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
