Eigen 创建求解稀疏矩阵 保存方便日后使用

来源:互联网 发布:c语言graphics 编辑:程序博客网 时间:2024/06/05 00:47
#include "Eigen/Sparse"  typedef Eigen::SparseMatrix<double> SparseMatrixType;    #include <vector>  using namespace std;  vec CDepthRegression::Get_LightDirection()  {      vec light;      int fn = m_mesh->faces.size();        typedef Eigen::Triplet<double> T;      std::vector<T> tripletList;//稀疏矩阵的元素      for (int i = 0; i < fn;i++)      {                    const TriMesh::Face &f = m_mesh->faces[i];          vec v0 = m_mesh->vertices[f[0]];          vec v1 = m_mesh->vertices[f[1]];          vec v2 = m_mesh->vertices[f[2]];            vec facenormal = (v2 - v1) CROSS(v0 - v2);          facenormal = normalize(facenormal);          for (int j = 0; j < 3;j++)          {              tripletList.push_back(T(i,j, facenormal[j]));          }                }      SparseMatrixType Ls(fn, 3);//矩阵的宽高      Ls.setFromTriplets(tripletList.begin(), tripletList.end());//从Triplet中构建稀疏矩阵        //最小二乘解超静定方程组      SparseMatrixType ls_transpose = Ls.transpose();      SparseMatrixType LsLs = ls_transpose* Ls;      Eigen::VectorXd RHSPos;//超静定方程组右边      RHSPos.resize(fn);      RHSPos.setZero();      for (int i = 0; i < fn;i++)      {          float I = m_mesh->face_color[i];          RHSPos[i] = I;      }        Eigen::SimplicialCholesky<SparseMatrixType>MatricesCholesky(LsLs);          Eigen::VectorXd xyzRHS = ls_transpose*RHSPos;      Eigen::Vector3d xyz = MatricesCholesky.solve(xyzRHS);          return vec(xyz[0],xyz[1], xyz[2]);    }