Libigl学习笔记——二、第二章——离散几何量和算子

最新推荐文章于 2024-05-01 21:12:48 发布
最新推荐文章于 2024-05-01 21:12:48 发布
阅读量449 收藏 1
点赞数 1
分类专栏： Libigl学习笔记 文章标签： 学习 笔记
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/weixin_45029825/article/details/132814562
版权

Libigl学习笔记——二、第二章——离散几何量和算子


本章说明了 libigl 可以在网格上计算的几个离散量，以及构造常用离散微分几何算子的 libigl 函数。它还介绍了查看器的基本绘图和着色程序。

第二章

本章中内容：

  • libigl 可以在网格上计算的几个离散量
  • 构造常用离散微分几何算子的 libigl 函数
  • 介绍了查看器的基本绘图和着色程序

2.1 网格的离散量

2.1.1 法线 Normals

  • 表面法线是渲染表面所需的基本量
  • 有多种方法可以在三角形网格上计算和存储法线
  • 例 201 演示如何使用 libigl 计算和可视化法线

2.1.2 每个面 Per-face

在网格的每个三角形上，法线被很好地定义为与三角形平面正交的矢量

  • 这些分段常数法线产生分段平坦渲染：表面看起来不光滑，并揭示了其潜在的离散化。

2.1.3 每个顶点 Per-vertex

  • 法线可以计算并存储在顶点上，并在三角形的内部插值以生成平滑渲染（Phong 着色）
  • 大多数计算每个顶点法线的技术都需要平均入射面法线
  • 主要区别在于它们的加权方案：
    • 均匀加权受到离散化选择的严重偏差
    • 而基于面积或基于角度的加权更宽容

基于面积的权重的典型半边样式计算具有以下结构：

N.setZero(V.rows(),3);
for(int i : vertices)
{
  for(face : incident_faces(i))
  {
    N.row(i) += face.area * face.normal;
  }
}
N.rowwise().normalize();
  • 在不使用半边数据结构的情况下循环事件面，从而构建每个顶点法线似乎效率低下。
  • 但是，每个顶点法线可能会将每个面垂直于其角顶点上的运行总和：
N.setZero(V.rows(),3);
for(int f = 0; f < F.rows();f++)
{
  for(int c = 0; c < 3;c++)
  {
    N.row(F(f,c)) += area(f) * face_normal.row(f);
  }
}
N.rowwise().normalize();

2.1.4 每个角 Per-corner

按角存储法线是支持平滑和锐利（例如折痕和拐角）渲染的一种高效便捷的方法

  • 这种格式对于OpenGL和.obj网格文件格式是通用的
  • 通常，网格设计者会调整此类法线，但折痕和拐角也可以自动计算。
  • Libigl 实现了一个简单的方案，该方案将角法线计算为入射到相应顶点上的面的法线的平均值，这些面的偏差不超过指定的二面角（例如 20°）
#include <igl/read_triangle_mesh.h>
#include <igl/opengl/glfw/Viewer.h>
#include <igl/per_vertex_normals.h>
#include <igl/per_face_normals.h>
#include <igl/per_corner_normals.h>
#include <iostream>

Eigen::MatrixXd V;
Eigen::MatrixXi F;

Eigen::MatrixXd N_vertices;
Eigen::MatrixXd N_faces;
Eigen::MatrixXd N_corners;


// This function is called every time a keyboard button is pressed
bool key_down(igl::opengl::glfw::Viewer& viewer, unsigned char key, int modifier)
{
  switch(key)
  {
    case '1':
      viewer.data().set_normals(N_faces);
      return true;
    case '2':
      viewer.data().set_normals(N_vertices);
      return true;
    case '3':
      viewer.data().set_normals(N_corners);
      return true;
    default: break;
  }
  return false;
}

int main(int argc, char *argv[])
{
  // Load a mesh in OFF format
  igl::read_triangle_mesh(
    argc>1?argv[1]: TUTORIAL_SHARED_PATH "/fandisk.off",V,F);

  // Compute per-face normals
  igl::per_face_normals(V,F,N_faces);

  // Compute per-vertex normals
  igl::per_vertex_normals(V,F,N_vertices);

  // Compute per-corner normals, |dihedral angle| > 20 degrees --> crease
  igl::per_corner_normals(V,F,20,N_corners);

  // Plot the mesh
  igl::opengl::glfw::Viewer viewer;
  viewer.callback_key_down = &key_down;
  viewer.data().show_lines = false;
  viewer.data().set_mesh(V, F);
  viewer.data().set_normals(N_faces);
  std::cout<<
    "Press '1' for per-face normals."<<std::endl<<
    "Press '2' for per-vertex normals."<<std::endl<<
    "Press '3' for per-corner normals."<<std::endl;
  viewer.launch();
}

在这里插入图片描述该 Normals 示例计算每个面（左）、每个顶点（中间）和每个角（右）的法线

2.2 高斯曲率 Gaussian Curvature

连续曲面上的高斯曲率定义为主曲率的乘积：

  • kG=k1k2

  • 作为一种内在度量，它取决于度量而不是表面的嵌入

  • 直观地说，高斯曲率告诉表面的局部球形或椭圆形程度 （ kG>0），表面的局部马鞍形或双曲线程度（ kG<0 ），或表面的局部圆柱形或抛物线（ kG=0 ）程度。

  • 在离散设置中，三角形网格上“离散高斯曲率”的一个定义是通过顶点的角度缺陷：
    在这里插入图片描述

  • N(i)是入射在三角形顶点上的位置， θij 是三角形 j 中顶点 i i 的角度 3

  • 像连续模拟一样，我们的离散高斯曲率揭示了域上的椭圆、双曲和抛物线顶点

  • 202完整展示：

#include <igl/gaussian_curvature.h>
#include <igl/massmatrix.h>
#include <igl/invert_diag.h>
#include <igl/readOFF.h>
#include <igl/opengl/glfw/Viewer.h>

int main(int argc, char *argv[])
{
  using namespace Eigen;
  using namespace std;
  MatrixXd V;
  MatrixXi F;
  igl::readOFF(TUTORIAL_SHARED_PATH "/bumpy.off",V,F);

  VectorXd K;
  // Compute integral of Gaussian curvature
  igl::gaussian_curvature(V,F,K);
  // Compute mass matrix
  SparseMatrix<double> M,Minv;
  igl::massmatrix(V,F,igl::MASSMATRIX_TYPE_DEFAULT,M);
  igl::invert_diag(M,Minv);
  // Divide by area to get integral average
  K = (Minv*K).eval();

  // Plot the mesh with pseudocolors
  igl::opengl::glfw::Viewer viewer;
  viewer.data().set_mesh(V, F);
  viewer.data().set_data(K);
  viewer.launch();
}

在这里插入图片描述该 GaussianCurvature 示例计算离散高斯曲率并以伪彩色将其可视化

2.3 曲率方向 Curvature Directions

  1. 曲面上某一点 (k1,k2) 处的两个主曲率测量曲面在不同方向上的弯曲程度。最大和最小（带符号）弯曲的方向称为主方向，并且始终是正交的。

  2. 平均曲率定义为主曲率的平均值：
    在这里插入图片描述

  3. 提取平均曲率的一种方法是检查应用于曲面位置的拉普拉斯-贝尔特拉米算子。结果是所谓的平均曲率法线：
    −Δx=Hn

  4. 使用余切拉普拉斯-贝尔特拉米算子 3 在 libigl 中的离散三角形网格上很容易计算这一点

#include <igl/cotmatrix.h>
#include <igl/massmatrix.h>
#include <igl/invert_diag.h>
...
MatrixXd HN;
SparseMatrix<double> L,M,Minv;
igl::cotmatrix(V,F,L);
igl::massmatrix(V,F,igl::MASSMATRIX_TYPE_VORONOI,M);
igl::invert_diag(M,Minv);
HN = -Minv*(L*V);
H = HN.rowwise().norm(); //up to sign
  1. 结合离散高斯曲率的角度缺陷定义，可以定义主曲率并使用最小二乘拟合来寻找方向 。
  2. 或者，确定主曲率的鲁棒方法是通过二次拟合。在每个顶点周围的邻域中，找到一个最佳拟合二次曲线，并在此二次曲线上分析计算主曲率值和方向（示例 203）。
  3. 示例203完整展示：
#include <igl/avg_edge_length.h>
#include <igl/cotmatrix.h>
#include <igl/invert_diag.h>
#include <igl/massmatrix.h>
#include <igl/parula.h>
#include <igl/per_corner_normals.h>
#include <igl/per_face_normals.h>
#include <igl/per_vertex_normals.h>
#include <igl/principal_curvature.h>
#include <igl/read_triangle_mesh.h>
#include <igl/opengl/glfw/Viewer.h>

Eigen::MatrixXd V;
Eigen::MatrixXi F;

int main(int argc, char *argv[])
{
  using namespace Eigen;
  std::string filename = TUTORIAL_SHARED_PATH "/fertility.off";
  if(argc>1)
  {
    filename = argv[1];
  }
  // Load a mesh in OFF format
  igl::read_triangle_mesh(filename, V, F);

  // Alternative discrete mean curvature
  MatrixXd HN;
  SparseMatrix<double> L,M,Minv;
  igl::cotmatrix(V,F,L);
  igl::massmatrix(V,F,igl::MASSMATRIX_TYPE_VORONOI,M);
  igl::invert_diag(M,Minv);
  // Laplace-Beltrami of position
  HN = -Minv*(L*V);
  // Extract magnitude as mean curvature
  VectorXd H = HN.rowwise().norm();

  // Compute curvature directions via quadric fitting
  MatrixXd PD1,PD2;
  VectorXd PV1,PV2;
  igl::principal_curvature(V,F,PD1,PD2,PV1,PV2);
  // mean curvature
  H = 0.5*(PV1+PV2);

  igl::opengl::glfw::Viewer viewer;
  viewer.data().set_mesh(V, F);

  viewer.data().set_data(H);

  // Average edge length for sizing
  const double avg = igl::avg_edge_length(V,F);

  // Draw a red segment parallel to the maximal curvature direction
  const RowVector3d red(0.8,0.2,0.2),blue(0.2,0.2,0.8);
  viewer.data().add_edges(V + PD1*avg, V - PD1*avg, red);

  // Draw a blue segment parallel to the minimal curvature direction
  viewer.data().add_edges(V + PD2*avg, V - PD2*avg, blue);

  // Hide wireframe
  viewer.data().show_lines = false;

  viewer.launch();
}

在这里插入图片描述该 CurvatureDirections 示例通过二次拟合计算主曲率，并使用交叉场可视化伪彩色和主方向的平均曲率

2.4 梯度 Gradient

  1. 表面上的标量函数可以离散化为分段线性函数，并在每个网格顶点处定义值：
    在这里插入图片描述

  2. 其中 ϕi 是由网格定义的分段线性帽子函数，因此对于每个三角形 ϕi 都是线性函数，该线性函数仅在顶点 i 处为一个，在其他角处为零。
    在这里插入图片描述帽子函数 ϕi 在顶点处为一个，在所有其他顶点 i 处为零，在入射三角形上为线性

  3. 因此，这种分段线性函数的梯度只是帽子函数梯度的总和：
    在这里插入图片描述

  4. 这表明梯度是 fi 值向量的线性函数。因为 在每个三角形中都是线性的，所以它们的梯度在每个 ϕi 三角形中都是恒定的。因此，我们的离散梯度运算符可以写成矩阵乘法，将顶点值转换为三角形值：
    在这里插入图片描述

  5. 其中 f 是 n×1 和 G 是一个 md×n 稀疏矩阵。该矩阵 G 可以通过几何方式推导，例如第 2 章 。Libigl grad 函数计算 G 三角形和四面体网格（示例 204）：

  6. 204案例展示：

#include <igl/avg_edge_length.h>
#include <igl/barycenter.h>
#include <igl/grad.h>
#include <igl/jet.h>
#include <igl/readDMAT.h>
#include <igl/readOFF.h>
#include <igl/opengl/glfw/Viewer.h>

#include <iostream>

int main(int argc, char *argv[])
{
  using namespace Eigen;
  using namespace std;
  MatrixXd V;
  MatrixXi F;

  // Load a mesh in OFF format
  igl::readOFF(TUTORIAL_SHARED_PATH "/cheburashka.off", V, F);

  // Read scalar function values from a file, U: #V by 1
  VectorXd U;
  igl::readDMAT(TUTORIAL_SHARED_PATH "/cheburashka-scalar.dmat",U);

  // Compute gradient operator: #F*3 by #V
  SparseMatrix<double> G;
  igl::grad(V,F,G);

  // Compute gradient of U
  MatrixXd GU = Map<const MatrixXd>((G*U).eval().data(),F.rows(),3);
  // Compute gradient magnitude
  const VectorXd GU_mag = GU.rowwise().norm();

  igl::opengl::glfw::Viewer viewer;
  viewer.data().set_mesh(V, F);

  viewer.data().set_data(U);

  // Average edge length divided by average gradient (for scaling)
  const double max_size = igl::avg_edge_length(V,F) / GU_mag.mean();
  // Draw a black segment in direction of gradient at face barycenters
  MatrixXd BC;
  igl::barycenter(V,F,BC);
  const RowVector3d black(0,0,0);
  viewer.data().add_edges(BC,BC+max_size*GU, black);

  // Hide wireframe
  viewer.data().show_lines = false;

  viewer.launch();
}

在这里插入图片描述该 Gradient 示例计算网格上输入函数的梯度并可视化向量场

2.5 拉普拉斯算子 Laplacian

  1. 离散拉普拉斯是一种必不可少的几何加工工具。拉普拉斯和拉普拉斯-贝尔特拉米算子存在许多解释和风味。

  2. 在开欧几里得空间中，拉普拉斯算子是梯度的通常发散（或者等价地，函数的拉普拉斯算子是其黑森的迹线）：
    在这里插入图片描述

  3. 拉普拉斯-贝尔特拉米算子将其推广到曲面。

  4. 在考虑三角形网格上的分段线性函数时，可以通过多种方式导出离散拉普拉斯函数。几何加工中最流行的是所谓的 余切拉普拉斯， L 它同时由FEM，DEC产生，并将散度定理应用于顶点一环。作为将顶点值转换为顶点值的线性算子，拉普拉斯矩阵是一个包含元素的 n×n 矩阵 L ：
    在这里插入图片描述

  5. 其中 N(i) ，与（相邻）顶点 i 相邻的顶点，以及与 αij,βij 边 ij 相对的角度。此公式导致用于构造的典型半边样式实现 L ：

for(int i : vertices)
{
  for(int j : one_ring(i))
  {
    for(int k : triangle_on_edge(i,j))
    {
      L(i,j) += cot(angle(i,j,k));
      L(i,i) -= cot(angle(i,j,k));
    }
  }
}
  1. 与以前类似，它似乎在没有半边数据结构的情况下循环单环。然而，事实并非如此，因为拉普拉斯算符可以通过将每个三角形的贡献相加来构建，这与狄利克雷能量的有限元离散化（平方梯度之和）非常吻合：
for(triangle t : triangles)
{
  for(edge i,j : t)
  {
    L(i,j) += cot(angle(i,j,k));
    L(j,i) += cot(angle(i,j,k));
    L(i,i) -= cot(angle(i,j,k));
    L(j,j) -= cot(angle(i,j,k));
  }
}
  1. Libigl 为三角形网格和四面体网格实现了离散的“余切”拉普拉斯量，使用快速几何规则构建两者，而不是“按书本”构建 FEM 构造，后者涉及许多（小）矩阵反演。
  2. 应用于网格顶点位置的算子相当于通过沿平均曲率法线方向流动表面来平滑（示例205）。请注意，这相当于最小化表面积。
  3. 205案例完整展示：
#include <igl/barycenter.h>
#include <igl/cotmatrix.h>
#include <igl/doublearea.h>
#include <igl/grad.h>
#include <igl/jet.h>
#include <igl/massmatrix.h>
#include <igl/per_vertex_normals.h>
#include <igl/readDMAT.h>
#include <igl/readOFF.h>
#include <igl/repdiag.h>
#include <igl/opengl/glfw/Viewer.h>

#include <iostream>

Eigen::MatrixXd V,U;
Eigen::MatrixXi F;
Eigen::SparseMatrix<double> L;
igl::opengl::glfw::Viewer viewer;

int main(int argc, char *argv[])
{
  using namespace Eigen;
  using namespace std;

  // Load a mesh in OFF format
  igl::readOFF(TUTORIAL_SHARED_PATH "/cow.off", V, F);

  // Compute Laplace-Beltrami operator: #V by #V
  igl::cotmatrix(V,F,L);

  // Alternative construction of same Laplacian
  SparseMatrix<double> G,K;
  // Gradient/Divergence
  igl::grad(V,F,G);
  // Diagonal per-triangle "mass matrix"
  VectorXd dblA;
  igl::doublearea(V,F,dblA);
  // Place areas along diagonal #dim times
  const auto & T = 1.*(dblA.replicate(3,1)*0.5).asDiagonal();
  // Laplacian K built as discrete divergence of gradient or equivalently
  // discrete Dirichelet energy Hessian
  K = -G.transpose() * T * G;
  cout<<"|K-L|: "<<(K-L).norm()<<endl;

  const auto &key_down = [](igl::opengl::glfw::Viewer &viewer,unsigned char key,int mod)->bool
  {
    switch(key)
    {
      case 'r':
      case 'R':
        U = V;
        break;
      case ' ':
      {
        // Recompute just mass matrix on each step
        SparseMatrix<double> M;
        igl::massmatrix(U,F,igl::MASSMATRIX_TYPE_BARYCENTRIC,M);
        // Solve (M-delta*L) U = M*U
        const auto & S = (M - 0.001*L);
        Eigen::SimplicialLLT<Eigen::SparseMatrix<double > > solver(S);
        assert(solver.info() == Eigen::Success);
        U = solver.solve(M*U).eval();
        // Compute centroid and subtract (also important for numerics)
        VectorXd dblA;
        igl::doublearea(U,F,dblA);
        double area = 0.5*dblA.sum();
        MatrixXd BC;
        igl::barycenter(U,F,BC);
        RowVector3d centroid(0,0,0);
        for(int i = 0;i<BC.rows();i++)
        {
          centroid += 0.5*dblA(i)/area*BC.row(i);
        }
        U.rowwise() -= centroid;
        // Normalize to unit surface area (important for numerics)
        U.array() /= sqrt(area);
        break;
      }
      default:
        return false;
    }
    // Send new positions, update normals, recenter
    viewer.data().set_vertices(U);
    viewer.data().compute_normals();
    viewer.core().align_camera_center(U,F);
    return true;
  };


  // Use original normals as pseudo-colors
  MatrixXd N;
  igl::per_vertex_normals(V,F,N);
  MatrixXd C = N.rowwise().normalized().array()*0.5+0.5;

  // Initialize smoothing with base mesh
  U = V;
  viewer.data().set_mesh(U, F);
  viewer.data().set_colors(C);
  viewer.callback_key_down = key_down;

  cout<<"Press [space] to smooth."<<endl;;
  cout<<"Press [r] to reset."<<endl;;
  return viewer.launch();
}

在这里插入图片描述该 Laplacian 示例使用余切拉普拉斯计算 2 共形平均曲率流

2.6 质量矩阵 Mass Matrix

  1. 质量矩阵是另一个 n×n 将顶点值转换为顶点值的矩阵 M 。
  • 从 FEM 的角度来看，它是内积的离散化：它考虑了每个顶点周围的面积。
  • 因此，通常是对角矩阵， M 使得 Mii 网格中顶点 i 周围的重心或沃罗诺伊区域 3 。
  • 该矩阵的逆矩阵也非常有用，因为它将积分量转换为逐点量，例如：
    在这里插入图片描述

  1. 通常，当遇到曲面上积分的平方量时，在顶点采样函数值时，质量矩阵将用作内积的离散化：
    在这里插入图片描述

  2. 另一种质量矩阵是将三角形向量值转换为三角形向量值的 md×md 矩阵 T 。该矩阵表示一个内积，说明与每个三角形相关的面积（即三角形真实面积）。

2.7 拉普拉斯的替代构造 Alternative Construction Of Laplacian

  1. 离散余切拉普拉斯的另一种构造是通过“平方”离散梯度算子。这可以通过应用格林恒等式（暂时忽略边界条件）得出：
    在这里插入图片描述

  2. 或者以矩阵形式，可立即转换为代码：
    在这里插入图片描述

  3. 所以我们有那个 L=GTTG .这也暗示了我们可以考虑 GT 为离散散度算子，因为拉普拉斯算子是梯度的散度。自然，是一个 n×md 稀疏矩阵， GT 它将存储在三角形面上的向量值转换为顶点处的标量散度值。

2.8 精确离散测地线距离 Exact Discrete Geodesic Distances

  1. 两点之间的离散测地线距离是限制在表面上的最短路径的长度。对于三角形网格，此类路径由一组段组成，这些段可以是网格的边或与三角形相交的。
  2. Libigl 包含一个由 Danil Kirsanov （https://code.google.com/archive/p/geodesic/） 开发的精确测地线算法的包装器，通过基于特征的 API 公开它。函数
igl::exact_geodesic(V,F,VS,FS,VT,FT,d);
  1. 计算每个顶点（VT）或面（FT）中最近的测地线距离，距输入网格 V，F 的源顶点 VS 或面 FS。
  • 输出写入向量 d，它首先列出 VT 中顶点的距离，然后列出 FT 中面的距离。
  • 例如，如果要计算从带有 id vid 的顶点到F 的所有顶点的距离，则可以使用：
Eigen::VectorXi VS,FS,VT,FT;
// The selected vertex is the source
VS.resize(1);
VS << vid;
// All vertices are the targets
VT.setLinSpaced(V.rows(),0,V.rows()-1);
Eigen::VectorXd d;
igl::exact_geodesic(V,F,VS,FS,VT,FT,d);
  1. 206案例完整展示
#include <igl/readOBJ.h>
#include <igl/opengl/glfw/Viewer.h>
#include <igl/exact_geodesic.h>
#include <igl/unproject_onto_mesh.h>
#include <igl/parula.h>
#include <igl/isolines_map.h>
#include <igl/PI.h>
#include <iostream>


int main(int argc, char *argv[])
{
  using namespace Eigen;
  using namespace std;
  Eigen::MatrixXd V;
  Eigen::MatrixXi F;
  igl::opengl::glfw::Viewer viewer;
  // Load a mesh in OFF format
  igl::readOBJ(TUTORIAL_SHARED_PATH "/armadillo.obj", V, F);

  const auto update_distance = [&](const int vid)
  {
    Eigen::VectorXi VS,FS,VT,FT;
    // The selected vertex is the source
    VS.resize(1);
    VS << vid;
    // All vertices are the targets
    VT.setLinSpaced(V.rows(),0,V.rows()-1);
    Eigen::VectorXd d;
    std::cout<<"Computing geodesic distance to vertex "<<vid<<"..."<<std::endl;
    igl::exact_geodesic(V,F,VS,FS,VT,FT,d);
    // Plot the mesh
    Eigen::MatrixXd CM;
    igl::parula(Eigen::VectorXd::LinSpaced(21,0,1).eval(),false,CM);	//21调整划分层数
    igl::isolines_map(Eigen::MatrixXd(CM),CM);
    viewer.data().set_colormap(CM);
    viewer.data().set_data(d);
  };

  // Plot a distance when a vertex is picked
  viewer.callback_mouse_down =
  [&](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))
    {
      int max;
      bc.maxCoeff(&max);
      int vid = F(fid,max);
      update_distance(vid);
      return true;
    }
    return false;
  };
  viewer.data().set_mesh(V,F);
  viewer.data().show_lines = false;

  cout << "Click on mesh to define new source.\n" << std::endl;
  update_distance(0);
  return viewer.launch();

允许交互式地选择源顶点并使用周期性颜色模式显示距离。
在这里插入图片描述在这里插入图片描述在这里插入图片描述在这里插入图片描述

2.9 多边形拉普拉斯算子 Polygon Laplactian

  1. 207案例完整展示
#include <igl/readOBJ.h>
#include <igl/cotmatrix.h>
#include <igl/massmatrix.h>
#include <igl/edges.h>
#include <igl/find.h>
#include <igl/min_quad_with_fixed.h>
#include <igl/polygons_to_triangles.h>
#include <igl/opengl/glfw/Viewer.h>
#include <iostream>

int main(int argc, char *argv[])
{
  using namespace Eigen;
  using namespace std;

  // Load a mesh in OBJ format
  Eigen::MatrixXd OV,V;
  Eigen::VectorXi I,C;
  igl::readOBJ(
    argc<=1?TUTORIAL_SHARED_PATH "/cylinder.obj" :argv[1],V,I,C);
  // center
  V.rowwise() -= V.colwise().mean();
  OV = V;
  // Convert polygon representation to triangles
  Eigen::MatrixXi F;
  Eigen::VectorXi J;
  igl::polygons_to_triangles(I,C,F,J);
  Eigen::SparseMatrix<double> pL,pM,pP;
  igl::cotmatrix(V,I,C,pL,pM,pP);
  Eigen::SparseMatrix<double> tL,tM;
  igl::cotmatrix(V,F,tL);
  igl::massmatrix(V,F,igl::MASSMATRIX_TYPE_DEFAULT,tM);
  const double bbd = (V.colwise().maxCoeff()- V.colwise().minCoeff()).norm();
  igl::opengl::glfw::Viewer vr;
  vr.data_list[0].set_mesh(V,F);
  vr.append_mesh();
  vr.selected_data_index = 0;
  vr.data_list[0].set_face_based(true);
  Eigen::MatrixXi E;
  igl::edges(I,C,E);
  bool show_edges = true;
  bool use_poly = true;
  Eigen::MatrixXd pHN;
  Eigen::MatrixXd tHN;
  const auto update = [&]()
  {
    // why does windows need this Eigen::VectorXd(...) ?
    pHN = (pL*V).array().colwise() / Eigen::VectorXd(pM.diagonal()).array();
    tHN = (tL*V).array().colwise() / Eigen::VectorXd(tM.diagonal()).array();
    pHN *= 1.0/pHN.rowwise().norm().maxCoeff();
    tHN *= 1.0/tHN.rowwise().norm().maxCoeff();
    const auto was_face_based  = vr.data_list[0].face_based;
    Eigen::MatrixXd QV(V.rows()*2,3);
    QV.topRows(V.rows()) = V;
    if(use_poly)
    {
      printf("using polygon Laplacian\n");
      QV.bottomRows(V.rows()) = V-pHN;
    }else
    {
      printf("using triangle Laplacian\n");
      QV.bottomRows(V.rows()) = V-tHN;
    }
    Eigen::MatrixXi QE(V.rows(),2);
    for(int i = 0;i<V.rows();i++){ QE(i,0)=i;QE(i,1)=i+V.rows();}
    vr.data_list[1].set_edges(QV,QE,Eigen::RowVector3d(1,1,1));

    if(use_poly)
    {
      vr.data_list[0].show_lines = false;
      if(show_edges)
      {
        vr.data_list[0].set_edges(V,E,Eigen::RowVector3d(0,0,0));
      }else
      {
        vr.data_list[0].clear_edges();
      }
    }else
    {
      vr.data_list[0].clear_edges();
      vr.data_list[0].show_lines = show_edges;
    }
    vr.data_list[0].set_face_based(was_face_based);
  };
  const double original_area = pM.diagonal().sum();
  const auto recompute_M = [&]()
  {
    Eigen::SparseMatrix<double> _1,_2;
    igl::cotmatrix(V,I,C,_1,pM,_2);
    igl::massmatrix(V,F,igl::MASSMATRIX_TYPE_DEFAULT,tM);
    V *= sqrt(original_area / pM.diagonal().sum());
  };
  const auto cmcf_step = [&]()
  {
    const Eigen::SparseMatrix<double> S =
      use_poly? ((pM) - 0.05*(pL)): ((tM) - 0.05*(tL));
    const Eigen::MatrixXd rhs = use_poly? pM*V : tM*V;
    Eigen::SimplicialLLT<Eigen::SparseMatrix<double > > solver(S);
    assert(solver.info() == Eigen::Success);
    V = solver.solve(rhs).eval();
    // recompute just mass matrices
    recompute_M();
    // center
    V.rowwise() -= V.colwise().mean();
    vr.data_list[0].set_vertices(V);
    vr.data_list[0].compute_normals();
    update();
  };
  vr.callback_key_pressed = [&](decltype(vr) &,unsigned int key, int mod)
  {
    switch(key)
    {
      case ' ': cmcf_step(); return true;
      case 'R': case 'r': V=OV;recompute_M();vr.data_list[0].set_vertices(V);vr.data_list[0].compute_normals(); update();return true;
      case 'P': case 'p': use_poly=!use_poly; update();return true;
      case 'L': case 'l': show_edges=!show_edges; update();return true;
    }
    return false;
  };
  update();
  vr.launch();
}
  • 多边形网格法线方向显示
    在这里插入图片描述
  • 多边形转化为三角形法线方向显示
    在这里插入图片描述
孟庆磊
关注
  • 1
    点赞
  • 踩
  • 1
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
专栏目录
拉普拉斯算子_Mesh is Art（9）：拉普拉斯算子
weixin_39581739的博客
12-16 3535
在计算机图形学中，曲面表达与处理一直是热门话题，它牵扯到几何建模、辅助设计等诸多重要应用。网格作为一种最流行的分段线性曲面表示方法，能够近似地表达原光滑曲面。本节我们将介绍网格的微分属性，它是根据某些定义在网格上的线性算子得出的，这些算子通常是网格的拉普拉斯算子的变体，提供了各种方便曲面表示与处理等应用的基础。相较于传统的被用于表达顶点位置的全局笛卡尔坐标，微分曲面表达承载着曲面局部形态...
散乱点求特定点高斯曲率和平均曲率
08-15
求散乱点中特定点的高斯曲率和平均曲率，使用点最近的五个点计算
matlab算法——计算三维散乱点云的曲率，包括主曲率，高斯曲率和平均曲率）
10-30
matlab算法计算三维散乱点云的曲率，包括主曲率，高斯曲率和平均曲率
openCV实战-系列教程6：图像金字塔与轮廓检测（图像金字塔/高斯金字塔/拉普拉斯金字塔/图像上采样/图像下采样/轮廓检测/轮廓特征与近似/模版匹配）、原理解析、源码解读
机器学习与软件工程
08-24 231
openCV实战-系列教程6：图像金字塔与轮廓检测（图像金字塔/高斯金字塔/拉普拉斯金字塔/图像上采样/图像下采样/轮廓检测/轮廓特征与近似/模版匹配）、原理解析、源码解读
几种网格平滑算法的实现
u013339596的专栏
02-13 1万+
网格平滑 网格平滑属于数字几何处理领域的问题，计算机图形学和计算机辅助设计中，用多边形网格可以表示复杂的三维实体。随着三维扫描和曲面重建技术的发展，得到这些实体表面的多边形网格表示已经不是难事，但所得到的表面往往包含含噪声。在形状设计领域，在散乱点拟合和光滑形伏、纹理映射等应用领域，都有对平滑曲面的极大需求。故产生了网格平滑这一个研究点。   平滑前
Libigl学习笔记——二、第八章——参考文献
weixin_45029825的博客
09-12 197
扫描曲面和体积的隐式建模，1994 年。阿米尔·瓦克斯曼、马塞尔·坎彭、奥尔加·迪亚曼蒂、丹尼尔·帕诺佐、大卫·博姆斯、克劳斯·希尔德布兰特、米蕾拉·本-陈。奥尔加·索金、亚龙·利普曼、丹尼尔·科恩-奥尔、马克·亚历克萨、克里斯蒂安·罗斯尔和汉斯-彼得·塞德尔。马蒂亚斯·埃克、托尼·德罗斯、汤姆·杜尚、休格斯·霍普、迈克尔·朗斯伯里、维尔纳·斯图茨尔。亚历克·雅各布森、伊利亚·巴兰、拉迪斯拉夫·卡万、约万·波波维奇和奥尔加·索尔金。加文·巴里尔、尼尔·迪克森、瑞安·施密特、大卫·莱文、亚历克·雅各布森。
【CVPR2019】论文完整列表一
热门推荐
TomRen
06-11 1万+
CVPR 2019 Paper list[1-1000] CVPR2019 完整列表二 论文题目与链接 Finding Task-Relevant Features for Few-Shot Learning by Category Traversal Edge-Labeling Graph Neural Network for Few-Shot Learning G...
三角网格表面高斯曲率的计算与可视化
fendahe的专栏
06-16 2613
三角网格表面高斯曲率的计算与可视化 好久没有写代码了，最近拿计算三角网格表面的高斯曲率练了练手，并实现了高斯曲率的可视化，复习了一点微分几何的知识。感觉有时候还是要自己把代码写出来，调试运行，结合试验结果，才能对相应的知识有更深的了解。 所谓曲面上某点的高斯曲率，即该点两个主曲率的乘积。把曲面上的顶点映射到单位球的球心,把法线的端点映射到球面上，即将曲面上的点与球面上的点建立了一种对应，叫做曲
Games102_lecture8几何建模与处理基础_离散微分几何，Utopia框架介绍
qq_41286360的博客
12-28 971
Games102_lecture8几何建模与处理基础_离散微分几何，Utopia框架介绍1 三角网格回顾2 三角网格编程初步3曲线曲面的微分几何1 曲线的几个概念2 曲面的几个概念4 离散微分几何1 连续估计2 离散估计5 极小曲面6 作业67 视频： 1 三角网格回顾 顶点，边 和面。 由边和顶点推出面， 但是由 边和面推不出顶点。 2 三角网格编程初步 原顶点，另外一个半边，所相邻的face，下一条边 3曲线曲面的微分几何 微分几何是研究一个曲面无穷小领域的导数，曲率的特征。 1 曲线的几个
Libigl学习笔记——二、第七章——其他
weixin_45029825的博客
09-12 1515
Libigl学习笔记——二、第七章——其他第七章7.1 网格统计 Mesh Statistics7.2 广义绕组数 Generalized Winding Number7.3 网格抽取 Mesh Decimation7.4 有符号距离 Signed Distances7.4.1 点位置 Point Location7.4.2 最近点 Closest Points7.4.3 符号距离 Signed Distance7.5 行进立方体 Marching Cubes7.6 小平面方向 Facet Orienta
libigl第七章-其他内容
weixin_42174695的博客
10-05 494
Libigl 包含各种各样的几何处理工具和函数，用于处理网格以及与其相关的线性代数：在本介绍性教程中讨论的内容太多了。我们在本章中提取了一些有趣的函数来强调。
离散微分几何和蛋白质复合物的结构研究
06-04
本文提出了一种新颖的三维方法进行蛋白质复合物的结构研究。 在该方法中，将使用一对二维三角锥（具有多个顶部顶点）的交点来描述蛋白质分子的表面。 作为蛋白质复合物的数学玩具模型，我们考虑n个单纯形（n = 2,3,4 ...）的闭合轨迹的复合物，其中蛋白质复合物的设计问题对应于汉密尔顿循环问题的扩展版本。 问题是找到n个单纯形的“一组”闭合轨迹，该轨迹填充由给定的一对n + 1维三角锥定义的n维区域。 在这里，我们使用三角形的离散微分几何（即2个单纯形）在n = 2的情况下给出扩展哈密顿循环问题的解决方案。
ComputationalGeometry:计算几何、离散微分几何、计算共形几何..
07-23
计算几何 欢迎 :waving_hand: 上帝永远几何化。 - 柏拉图 欢迎来到这个回购。 随着几何深度学习和几何处理的发展，几何在建筑、计算设计和机器人制造领域发挥着越来越重要的作用。 虽然我能够时不时地编写几何算法，但我的目标是对这个主题有一个全面而全面的了解。 因此，我觉得有必要研究不同学科的计算几何。 结构 :puzzle_piece: 针对不同的资料来源，本次repo整理如下： 书 :books: 必须读的书 演讲 :school: 流行的在线讲座 图书馆 :laptop: 流行的几何处理库 纸 :page_with_curl: 必读的论文 作坊 :video_camera: 工期短 推荐材料 :thumbs_up: 我推荐的学习材料如下： 名称 类型 关联 离散微分几何，Keenan Crane 在线讲座 数字几何处理，付晓明 在线讲座（中文） 计算共形几何 在线讲座 计算几何：算法和应用 书 学习进度 离散微分几何： 第01章概述 第02章组合曲面 第03章微分几何导论 第04章外微积分导论 第05章离散曲面
人人学IoT———第二章学习笔记
01-06
第二个挑战：标准不同一、集成困难 终端传感器厂家多，通信协议众多、标准不统一、集成困难等 第三个挑战：网络连接复杂 物联网环境复杂、网络类型多、厂家开发终端需要适配很多协议、厂家进入物联网领域门槛较高等 ...
Solidworks学习笔记——随形变化.docx
07-13
Solidworks学习笔记——随形变化.docxSolidworks学习笔记——随形变化.docxSolidworks学习笔记——随形变化.docxSolidworks学习笔记——随形变化.docxSolidworks学习笔记——随形变化.docxSolidworks学习笔记——...
c2.zip c语言学习第二章
07-05
c语言学习第二章 数据类型，变量，标识符规范等总结 和源码
Android学习笔记——Menu介绍(二)
01-04
这次将继续上一篇文章没有讲完的Menu的学习，上下文菜单(Context menu)和弹出菜单(Popup menu)。 上下文菜单 上下文菜单提供对UI界面上的特定项或上下文框架的操作，就如同Windows中右键菜单一样。 在Android中，有...
JBPM学习笔记——流程设计与控制
03-03
火龙果软件工程技术中心 相关资料：《jBPM学习笔记(V3.2环境部署)》《jBPM学习笔记(框架设计简介)》背景本片文章，我们将从业务流程的设计开始，通过带领大家完成一个完整工作流的程序设计，来学习jPDL的使用。...
限流的学习
最新发布
weixin_43675252的博客
05-01 314
限流算法：滑动窗口算法滑动日志算法漏桶算法令牌桶算法。
cda备考学习学习笔记——基础知识篇(二)
09-18
cda备考学习学习笔记——基础知识篇(二)主要涉及了计算机科学与技术领域的基本概念和知识。 首先，它介绍了计算机网络的基础知识。网络是将多台计算机通过通信链路连接起来，使它们能够相互通信和共享资源的系统。...

“相关推荐”对你有帮助么？

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值