C++ SparseMatrix::reserve方法代码示例

IGL_INLINE void igl::adjacency_matrix(
  const Eigen::PlainObjectBase<DerivedF> & F, 
  Eigen::SparseMatrix<T>& A)
{
  using namespace std;
  using namespace Eigen;
  typedef typename DerivedF::Scalar Index;

  typedef Triplet<T> IJV;
  vector<IJV > ijv;
  ijv.reserve(F.size()*2);
  // Loop over faces
  for(int i = 0;i<F.rows();i++)
  {
    // Loop over this face
    for(int j = 0;j<F.cols();j++)
    {
      // Get indices of edge: s --> d
      Index s = F(i,j);
      Index d = F(i,(j+1)%F.cols());
      ijv.push_back(IJV(s,d,1));
      ijv.push_back(IJV(d,s,1));
    }
  }

  const Index n = F.maxCoeff()+1;
  A.resize(n,n);
  switch(F.cols())
  {
    case 3:
      A.reserve(6*(F.maxCoeff()+1));
      break;
    case 4:
      A.reserve(26*(F.maxCoeff()+1));
      break;
  }
  A.setFromTriplets(ijv.begin(),ijv.end());

  // Force all non-zeros to be one

  // Iterate over outside
  for(int k=0; k<A.outerSize(); ++k)
  {
    // Iterate over inside
    for(typename Eigen::SparseMatrix<T>::InnerIterator it (A,k); it; ++it)
    {
      assert(it.value() != 0);
      A.coeffRef(it.row(),it.col()) = 1;
    }
  }
}

void LaplacianOperator::computeLaplacianOperator( Eigen::SparseMatrix<double>& laplacianOperator )
{
	laplacianOperator.resize(mMeshVertexCount,mMeshVertexCount);
	laplacianOperator.reserve(Eigen::VectorXi::Constant(mMeshVertexCount,10));
	for (int i = 0; i < mMeshVertexCount; i++)
	{
		/* 如果第i个点没有邻接点，即它是一个孤立的点，那么它的laplacian坐标为0 */
		if( adjacentMatrix.innerVector(i).nonZeros() == 0)
		{
			laplacianOperator.insert(i,i) = 0;
			continue;
		}
		laplacianOperator.insert(i,i) = 1;
#ifdef MY_DEBUG
		int adjCount = 0;
#endif		
		for (Eigen::SparseMatrix<double>::InnerIterator it(adjacentMatrix,i); it; it++)
		{
			if(i != it.row())
			{
				laplacianOperator.insert(i,it.row()) = -1/degreeMatrix(i);				
#ifdef MY_DEBUG
			adjCount++;
			if(adjCount >= 10)
				printf("InnerVector size should expand! CurrentMax:%d.\n",adjCount);
#endif
			}
		}
	}
}

IGL_INLINE void igl::speye(const int m, const int n, Eigen::SparseMatrix<T> & I)
{
  // size of diagonal
  int d = (m<n?m:n);
  I = Eigen::SparseMatrix<T>(m,n);
  I.reserve(d);
  for(int i = 0;i<d;i++)
  {
    I.insert(i,i) = 1.0;
  }
  I.finalize();
}

inline
void
identity_operator(
    Eigen::SparseMatrix<double>& result,
    int size)
{
  result.resize(size, size);
  result.reserve(size);
  std::vector<double_triplet_t> matrix_coeffs;
  matrix_coeffs.reserve(size); // #diag
  // diag
  for (int i = 0; i < size; ++i) matrix_coeffs.push_back(double_triplet_t(i,i, 1.));
  result.setFromTriplets(matrix_coeffs.begin(), matrix_coeffs.end());
}

IGL_INLINE void igl::cotmatrix(
  const Eigen::MatrixBase<DerivedV> & V, 
  const Eigen::MatrixBase<DerivedF> & F, 
  Eigen::SparseMatrix<Scalar>& L)
{
  using namespace Eigen;
  using namespace std;

  L.resize(V.rows(),V.rows());
  Matrix<int,Dynamic,2> edges;
  int simplex_size = F.cols();
  // 3 for triangles, 4 for tets
  assert(simplex_size == 3 || simplex_size == 4);
  if(simplex_size == 3)
  {
    // This is important! it could decrease the comptuation time by a factor of 2
    // Laplacian for a closed 2d manifold mesh will have on average 7 entries per
    // row
    L.reserve(10*V.rows());
    edges.resize(3,2);
    edges << 
      1,2,
      2,0,
      0,1;
  }else if(simplex_size == 4)
  {
    L.reserve(17*V.rows());
    edges.resize(6,2);
    edges << 
      1,2,
      2,0,
      0,1,
      3,0,
      3,1,
      3,2;
  }else
  {
    return;
  }
  // Gather cotangents
  Matrix<Scalar,Dynamic,Dynamic> C;
  cotmatrix_entries(V,F,C);
  
  vector<Triplet<Scalar> > IJV;
  IJV.reserve(F.rows()*edges.rows()*4);
  // Loop over triangles
  for(int i = 0; i < F.rows(); i++)
  {
    // loop over edges of element
    for(int e = 0;e<edges.rows();e++)
    {
      int source = F(i,edges(e,0));
      int dest = F(i,edges(e,1));
      IJV.push_back(Triplet<Scalar>(source,dest,C(i,e)));
      IJV.push_back(Triplet<Scalar>(dest,source,C(i,e)));
      IJV.push_back(Triplet<Scalar>(source,source,-C(i,e)));
      IJV.push_back(Triplet<Scalar>(dest,dest,-C(i,e)));
    }
  }
  L.setFromTriplets(IJV.begin(),IJV.end());
}

IGL_INLINE void igl::arap_linear_block_spokes(
  const MatV & V,
  const MatF & F,
  const int d,
  Eigen::SparseMatrix<Scalar> & Kd)
{
  using namespace std;
  using namespace Eigen;
  // simplex size (3: triangles, 4: tetrahedra)
  int simplex_size = F.cols();
  // Number of elements
  int m = F.rows();
  // Temporary output
  Matrix<int,Dynamic,2> edges;
  Kd.resize(V.rows(), V.rows());
  vector<Triplet<Scalar> > Kd_IJV;
  if(simplex_size == 3)
  {
    // triangles
    Kd.reserve(7*V.rows());
    Kd_IJV.reserve(7*V.rows());
    edges.resize(3,2);
    edges << 
      1,2,
      2,0,
      0,1;
  }else if(simplex_size == 4)
  {
    // tets
    Kd.reserve(17*V.rows());
    Kd_IJV.reserve(17*V.rows());
    edges.resize(6,2);
    edges << 
      1,2,
      2,0,
      0,1,
      3,0,
      3,1,
      3,2;
  }
  // gather cotangent weights
  Matrix<Scalar,Dynamic,Dynamic> C;
  cotmatrix_entries(V,F,C);
  // should have weights for each edge
  assert(C.cols() == edges.rows());
  // loop over elements
  for(int i = 0;i<m;i++)
  {
    // loop over edges of element
    for(int e = 0;e<edges.rows();e++)
    {
      int source = F(i,edges(e,0));
      int dest = F(i,edges(e,1));
      double v = 0.5*C(i,e)*(V(source,d)-V(dest,d));
      Kd_IJV.push_back(Triplet<Scalar>(source,dest,v));
      Kd_IJV.push_back(Triplet<Scalar>(dest,source,-v));
      Kd_IJV.push_back(Triplet<Scalar>(source,source,v));
      Kd_IJV.push_back(Triplet<Scalar>(dest,dest,-v));
    }
  }
  Kd.setFromTriplets(Kd_IJV.begin(),Kd_IJV.end());
  Kd.makeCompressed();
}

//.........这里部分代码省略.........

  const auto & apply_sort = [] (
     const MatrixX4I & T, 
     const MatrixX4I & sJ, 
     MatrixX4I & sT)
  {
    sT.resize(T.rows(),4);
    for(size_t t = 0;t<(size_t)T.rows();t++)
    {
      for(size_t c = 0;c<4;c++)
      {
        sT(t,c) = T(t,sJ(t,c));
      }
    }
  };

  const auto & one_below = [&V,&apply_sort](
    const MatrixX4I & T,
    const MatrixX4S & IT,
    MatrixX3I & G,
    SparseMatrix<BCType> & BC)
  {
    // Number of tets
    const size_t m = T.rows();
    MatrixX4S sIT;
    MatrixX4I sJ;
    sort(IT,2,true,sIT,sJ);
    MatrixX4I sT;
    apply_sort(T,sJ,sT);
    MatrixX3S lambda = 
      sIT.rightCols(3).array() /
      (sIT.rightCols(3).colwise()-sIT.col(0)).array();
    vector<Triplet<BCType> > IJV;
    IJV.reserve(m*3*2);
    for(size_t c = 0;c<3;c++)
    {
      for(size_t t = 0;t<(size_t)m;t++)
      {
        IJV.push_back(Triplet<BCType>(c*m+t,  sT(t,0),  lambda(t,c)));
        IJV.push_back(Triplet<BCType>(c*m+t,sT(t,c+1),1-lambda(t,c)));
      }
    }
    BC.resize(m*3,V.rows());
    BC.reserve(m*3*2);
    BC.setFromTriplets(IJV.begin(),IJV.end());
    G.resize(m,3);
    for(size_t c = 0;c<3;c++)
    {
      G.col(c).setLinSpaced(m,0+c*m,(m-1)+c*m);
    }
  };

  const auto & two_below = [&V,&apply_sort](
    const MatrixX4I & T,
    const MatrixX4S & IT,
    MatrixX3I & G,
    SparseMatrix<BCType> & BC)
  {
    // Number of tets
    const size_t m = T.rows();
    MatrixX4S sIT;
    MatrixX4I sJ;
    sort(IT,2,true,sIT,sJ);
    MatrixX4I sT;
    apply_sort(T,sJ,sT);
    MatrixX2S lambda =