C++ PlainObjectBase::block方法代码示例

IGL_INLINE void igl::fit_rotations_planar(
  const Eigen::PlainObjectBase<DerivedS> & S,
        Eigen::PlainObjectBase<DerivedD> & R)
{ 
  using namespace std;
  const int dim = S.cols();
  const int nr = S.rows()/dim;
  //assert(dim == 2 && "_planar input should be 2D");
  assert(nr * dim == S.rows());

  // resize output
  R.resize(dim,dim*nr); // hopefully no op (should be already allocated)

  Eigen::Matrix<typename DerivedS::Scalar,2,2> si;
  // loop over number of rotations we're computing
  for(int r = 0;r<nr;r++)
  {
    // build this covariance matrix
    for(int i = 0;i<2;i++)
    {
      for(int j = 0;j<2;j++)
      {
        si(i,j) = S(i*nr+r,j);
      }
    }
    typedef Eigen::Matrix<typename DerivedD::Scalar,2,2> Mat2;
    typedef Eigen::Matrix<typename DerivedD::Scalar,2,1> Vec2;
    Mat2 ri,ti,ui,vi;
    Vec2 _;
    igl::polar_svd(si,ri,ti,ui,_,vi);
#ifndef FIT_ROTATIONS_ALLOW_FLIPS
    // Check for reflection
    if(ri.determinant() < 0)
    {
      vi.col(1) *= -1.;
      ri = ui * vi.transpose();
    }
    assert(ri.determinant() >= 0);
#endif  

    // Not sure why polar_dec computes transpose...
    R.block(0,r*dim,dim,dim).setIdentity();
    R.block(0,r*dim,2,2) = ri.transpose();
  }
}

IGL_INLINE void igl::false_barycentric_subdivision(
    const Eigen::PlainObjectBase<Scalar> & V,
    const Eigen::PlainObjectBase<Index> & F,
    Eigen::PlainObjectBase<Scalar> & VD,
    Eigen::PlainObjectBase<Index> & FD)
{
  using namespace Eigen;
  // Compute face barycenter
  Eigen::MatrixXd BC;
  igl::barycenter(V,F,BC);

  // Add the barycenters to the vertices
  VD.resize(V.rows()+F.rows(),3);
  VD.block(0,0,V.rows(),3) = V;
  VD.block(V.rows(),0,F.rows(),3) = BC;

  // Each face is split four ways
  FD.resize(F.rows()*3,3);

  for (unsigned i=0; i<F.rows(); ++i)
  {
    int i0 = F(i,0);
    int i1 = F(i,1);
    int i2 = F(i,2);
    int i3 = V.rows() + i;

    Vector3i F0,F1,F2;
    F0 << i0,i1,i3;
    F1 << i1,i2,i3;
    F2 << i2,i0,i3;

    FD.row(i*3 + 0) = F0;
    FD.row(i*3 + 1) = F1;
    FD.row(i*3 + 2) = F2;
  }


}

IGL_INLINE void igl::fit_rotations(
  const Eigen::PlainObjectBase<DerivedS> & S,
  const bool single_precision,
  Eigen::PlainObjectBase<DerivedD> & R)
{
  using namespace std;
  const int dim = S.cols();
  const int nr = S.rows()/dim;
  assert(nr * dim == S.rows());
  assert(dim == 3);

  // resize output
  R.resize(dim,dim*nr); // hopefully no op (should be already allocated)

  //std::cout<<"S=["<<std::endl<<S<<std::endl<<"];"<<std::endl;
  //MatrixXd si(dim,dim);
  Eigen::Matrix<typename DerivedS::Scalar,3,3> si;// = Eigen::Matrix3d::Identity();
  // loop over number of rotations we're computing
  for(int r = 0;r<nr;r++)
  {
    // build this covariance matrix
    for(int i = 0;i<dim;i++)
    {
      for(int j = 0;j<dim;j++)
      {
        si(i,j) = S(i*nr+r,j);
      }
    }
    typedef Eigen::Matrix<typename DerivedD::Scalar,3,3> Mat3;
    typedef Eigen::Matrix<typename DerivedD::Scalar,3,1> Vec3;
    Mat3 ri;
    if(single_precision)
    {
      polar_svd3x3(si, ri);
    }else
    {
      Mat3 ti,ui,vi;
      Vec3 _;
      igl::polar_svd(si,ri,ti,ui,_,vi);
    }
    assert(ri.determinant() >= 0);
    R.block(0,r*dim,dim,dim) = ri.block(0,0,dim,dim).transpose();
    //cout<<matlab_format(si,C_STR("si_"<<r))<<endl;
    //cout<<matlab_format(ri.transpose().eval(),C_STR("ri_"<<r))<<endl;
  }
}

IGL_INLINE double igl::polyvector_field_poisson_reconstruction(
  const Eigen::PlainObjectBase<DerivedV> &Vcut,
  const Eigen::PlainObjectBase<DerivedF> &Fcut,
  const Eigen::PlainObjectBase<DerivedS> &sol3D_combed,
  Eigen::PlainObjectBase<DerivedSF> &scalars,
  Eigen::PlainObjectBase<DerivedS> &sol3D_recon,
  Eigen::PlainObjectBase<DerivedE> &max_error )
  {
    Eigen::SparseMatrix<typename DerivedV::Scalar> gradMatrix;
    igl::grad(Vcut, Fcut, gradMatrix);

    Eigen::VectorXd FAreas;
    igl::doublearea(Vcut, Fcut, FAreas);
    FAreas = FAreas.array() * .5;

    int nf = FAreas.rows();
    Eigen::SparseMatrix<typename DerivedV::Scalar> M,M1;
    Eigen::VectorXi II = igl::colon<int>(0, nf-1);

    igl::sparse(II, II, FAreas, M1);
    igl::repdiag(M1, 3, M) ;

    int half_degree = sol3D_combed.cols()/3;

    sol3D_recon.setZero(sol3D_combed.rows(),sol3D_combed.cols());

    int numF = Fcut.rows();
    scalars.setZero(Vcut.rows(),half_degree);

    Eigen::SparseMatrix<typename DerivedV::Scalar> Q = gradMatrix.transpose()* M *gradMatrix;

    //fix one point at Ik=fix, value at fixed xk=0
    int fix = 0;
    Eigen::VectorXi Ik(1);Ik<<fix;
    Eigen::VectorXd xk(1);xk<<0;

    //unknown indices
    Eigen::VectorXi Iu(Vcut.rows()-1,1);
    Iu<<igl::colon<int>(0, fix-1),  igl::colon<int>(fix+1,Vcut.rows()-1);

    Eigen::SparseMatrix<typename DerivedV::Scalar> Quu, Quk;
    igl::slice(Q, Iu, Iu, Quu);
    igl::slice(Q, Iu, Ik, Quk);
    Eigen::SimplicialLDLT<Eigen::SparseMatrix<typename DerivedV::Scalar> > solver;
    solver.compute(Quu);


    Eigen::VectorXd vec; vec.setZero(3*numF,1);
    for (int i =0; i<half_degree; ++i)
    {
      vec<<sol3D_combed.col(i*3+0),sol3D_combed.col(i*3+1),sol3D_combed.col(i*3+2);
      Eigen::VectorXd b = gradMatrix.transpose()* M * vec;
      Eigen::VectorXd bu = igl::slice(b, Iu);

      Eigen::VectorXd rhs = bu-Quk*xk;
      Eigen::MatrixXd yu = solver.solve(rhs);

      Eigen::VectorXi index = i*Eigen::VectorXi::Ones(Iu.rows(),1);
      igl::slice_into(yu, Iu, index, scalars);scalars(Ik[0],i)=xk[0];
    }

    //    evaluate gradient of found scalar function
    for (int i =0; i<half_degree; ++i)
    {
      Eigen::VectorXd vec_poisson = gradMatrix*scalars.col(i);
      sol3D_recon.col(i*3+0) = vec_poisson.segment(0*numF, numF);
      sol3D_recon.col(i*3+1) = vec_poisson.segment(1*numF, numF);
      sol3D_recon.col(i*3+2) = vec_poisson.segment(2*numF, numF);
    }

    max_error.setZero(numF,1);
    for (int i =0; i<half_degree; ++i)
    {
      Eigen::VectorXd diff = (sol3D_recon.block(0, i*3, numF, 3)-sol3D_combed.block(0, i*3, numF, 3)).rowwise().norm();
      diff = diff.array() / sol3D_combed.block(0, i*3, numF, 3).rowwise().norm().array();
      max_error = max_error.cwiseMax(diff.cast<typename DerivedE::Scalar>());
    }

    return max_error.mean();
  }

IGL_INLINE void igl::all_edges(
  const Eigen::PlainObjectBase<DerivedF> & F,
  Eigen::PlainObjectBase<DerivedE> & E)
{
  E.resize(F.rows()*F.cols(),F.cols()-1);
  switch(F.cols())
  {
    case 4:
      E.block(0*F.rows(),0,F.rows(),1) = F.col(1);
      E.block(0*F.rows(),1,F.rows(),1) = F.col(3);
      E.block(0*F.rows(),2,F.rows(),1) = F.col(2);

      E.block(1*F.rows(),0,F.rows(),1) = F.col(0);
      E.block(1*F.rows(),1,F.rows(),1) = F.col(2);
      E.block(1*F.rows(),2,F.rows(),1) = F.col(3);

      E.block(2*F.rows(),0,F.rows(),1) = F.col(0);
      E.block(2*F.rows(),1,F.rows(),1) = F.col(3);
      E.block(2*F.rows(),2,F.rows(),1) = F.col(1);

      E.block(3*F.rows(),0,F.rows(),1) = F.col(0);
      E.block(3*F.rows(),1,F.rows(),1) = F.col(1);
      E.block(3*F.rows(),2,F.rows(),1) = F.col(2);
      return;
    case 3:
      E.block(0*F.rows(),0,F.rows(),1) = F.col(1);
      E.block(0*F.rows(),1,F.rows(),1) = F.col(2);
      E.block(1*F.rows(),0,F.rows(),1) = F.col(2);
      E.block(1*F.rows(),1,F.rows(),1) = F.col(0);
      E.block(2*F.rows(),0,F.rows(),1) = F.col(0);
      E.block(2*F.rows(),1,F.rows(),1) = F.col(1);
      return;
  }
}

IGL_INLINE void igl::upsample(
  const Eigen::PlainObjectBase<DerivedV>& V,
  const Eigen::PlainObjectBase<DerivedF>& F,
  Eigen::PlainObjectBase<DerivedNV>& NV,
  Eigen::PlainObjectBase<DerivedNF>& NF)
{
  // Use "in place" wrapper instead
  assert(&V != &NV);
  assert(&F != &NF);
  using namespace igl;
  using namespace std;
  using namespace Eigen;

  ////typedef Eigen::PlainObjectBase<DerivedF> MatF;
  ////typedef Eigen::PlainObjectBase<DerivedNF> MatNF;
  ////typedef Eigen::PlainObjectBase<DerivedNV> MatNV;
  ////MatF FF, FFi;
  Eigen::MatrixXi FF,FFi;
  tt(V,F,FF,FFi);

  // TODO: Cache optimization missing from here, it is a mess
  
  // Compute the number and positions of the vertices to insert (on edges)
  Eigen::MatrixXi NI = Eigen::MatrixXi::Constant(FF.rows(),FF.cols(),-1);
  int counter = 0;
  
  for(int i=0;i<FF.rows();++i)
  {
    for(int j=0;j<3;++j)
    {
      if(NI(i,j) == -1)
      {
        NI(i,j) = counter;
        if (FF(i,j) != -1) // If it is not a border
          NI(FF(i,j),FFi(i,j)) = counter;
        ++counter;
      }
    }
  }

  int n_odd = V.rows();
  int n_even = counter;

  // Not sure what this is
  //Eigen::DynamicSparseMatrix<double> SUBD(V.rows()+n_even,V.rows());
  //SUBD.reserve(15 * (V.rows()+n_even));
  
  // Preallocate NV and NF
  NV.resize(V.rows()+n_even,V.cols());
  NF.resize(F.rows()*4,3);
  
  // Fill the odd vertices position
  NV.block(0,0,V.rows(),V.cols()) = V;

  // Fill the even vertices position
  for(int i=0;i<FF.rows();++i)
  {
    for(int j=0;j<3;++j)
    {
      NV.row(NI(i,j) + n_odd) = 0.5 * V.row(F(i,j)) + 0.5 * V.row(F(i,(j+1)%3));
    }
  }

  // Build the new topology (Every face is replaced by four)
  for(int i=0; i<F.rows();++i)
  {
    VectorXi VI(6);
    VI << F(i,0), F(i,1), F(i,2), NI(i,0) + n_odd, NI(i,1) + n_odd, NI(i,2) + n_odd;
    
    VectorXi f0(3), f1(3), f2(3), f3(3);
    f0 << VI(0), VI(3), VI(5);
    f1 << VI(1), VI(4), VI(3);
    f2 << VI(3), VI(4), VI(5);
    f3 << VI(4), VI(2), VI(5);
    
    NF.row((i*4)+0) = f0;
    NF.row((i*4)+1) = f1;
    NF.row((i*4)+2) = f2;
    NF.row((i*4)+3) = f3;
  }
  
}

//.........这里部分代码省略.........
  {
    const auto & bounding_box = [](
      const Eigen::PlainObjectBase<DerivedV> & V,
      const MatrixXG & F)->
        DerivedV
    {
      DerivedV BB(2,3);
      BB<<
         1e26,1e26,1e26,
        -1e26,-1e26,-1e26;
      const size_t m = F.rows();
      for(size_t f = 0;f<m;f++)
      {
        for(size_t c = 0;c<3;c++)
        {
          const auto & vfc = V.row(F(f,c)).eval();
          BB(0,0) = std::min(BB(0,0), vfc(0,0));
          BB(0,1) = std::min(BB(0,1), vfc(0,1));
          BB(0,2) = std::min(BB(0,2), vfc(0,2));
          BB(1,0) = std::max(BB(1,0), vfc(0,0));
          BB(1,1) = std::max(BB(1,1), vfc(0,1));
          BB(1,2) = std::max(BB(1,2), vfc(0,2));
        }
      }
      return BB;
    };
    // A lot of the time we're dealing with unrelated, distant components: cull
    // them.
    DerivedV ABB = bounding_box(V,A);
    DerivedV BBB = bounding_box(V,B);
    if( (BBB.row(0)-ABB.row(1)).maxCoeff()>0  ||
        (ABB.row(0)-BBB.row(1)).maxCoeff()>0 )
    {
      // bounding boxes do not overlap
      return false;
    } else {
      return true;
    }
  };

  // Reject components which are completely inside other components
  vector<bool> keep(ncc,true);
  size_t nG = 0;
  // This is O( ncc * ncc * m)
  for(size_t id = 0;id<ncc;id++)
  {
    if (!keep[id]) continue;
    std::vector<size_t> unresolved;
    for(size_t oid = 0;oid<ncc;oid++)
    {
      if(id == oid || !keep[oid])
      {
        continue;
      }
      if (has_overlapping_bbox(V, vG[id], vG[oid])) {
          unresolved.push_back(oid);
      }
    }
    const size_t num_unresolved_components = unresolved.size();
    DerivedV query_points(num_unresolved_components, 3);
    for (size_t i=0; i<num_unresolved_components; i++) {
        const size_t oid = unresolved[i];
        DerivedF f = vG[oid].row(0);
        query_points(i,0) = (V(f(0,0), 0) + V(f(0,1), 0) + V(f(0,2), 0))/3.0;
        query_points(i,1) = (V(f(0,0), 1) + V(f(0,1), 1) + V(f(0,2), 1))/3.0;
        query_points(i,2) = (V(f(0,0), 2) + V(f(0,1), 2) + V(f(0,2), 2))/3.0;
    }
    Eigen::VectorXi inside;
    igl::copyleft::cgal::points_inside_component(V, vG[id], query_points, inside);
    assert((size_t)inside.size() == num_unresolved_components);
    for (size_t i=0; i<num_unresolved_components; i++) {
        if (inside(i, 0)) {
            const size_t oid = unresolved[i];
            keep[oid] = false;
        }
    }
  }
  for (size_t id = 0; id<ncc; id++) {
      if (keep[id]) {
          nG += vJ[id].rows();
      }
  }

  // collect G and J across components
  G.resize(nG,3);
  J.resize(nG,1);
  {
    size_t off = 0;
    for(Index id = 0;id<(Index)ncc;id++)
    {
      if(keep[id])
      {
        assert(vG[id].rows() == vJ[id].rows());
        G.block(off,0,vG[id].rows(),vG[id].cols()) = vG[id];
        J.block(off,0,vJ[id].rows(),vJ[id].cols()) = vJ[id];
        off += vG[id].rows();
      }
    }
  }
}

//.........这里部分代码省略.........
    {
      MatrixXV BB(2,3);
      BB<<
         1e26,1e26,1e26,
        -1e26,-1e26,-1e26;
      const size_t m = F.rows();
      for(size_t f = 0;f<m;f++)
      {
        for(size_t c = 0;c<3;c++)
        {
          const auto & vfc = V.row(F(f,c));
          BB.row(0) = BB.row(0).array().min(vfc.array()).eval();
          BB.row(1) = BB.row(1).array().max(vfc.array()).eval();
        }
      }
      return BB;
    };
    // A lot of the time we're dealing with unrelated, distant components: cull
    // them.
    MatrixXV ABB = bounding_box(V,A);
    MatrixXV BBB = bounding_box(V,B);
    if( (BBB.row(0)-ABB.row(1)).maxCoeff()>0  ||
        (ABB.row(0)-BBB.row(1)).maxCoeff()>0 )
    {
      // bounding boxes do not overlap
      return false;
    }
    ////////////////////////////////////////////////////////////////////////
    // POTENTIAL ROBUSTNESS WEAK AREA
    ////////////////////////////////////////////////////////////////////////
    //
    // q could be so close (<~1e-16) to B that the winding number is not a robust way to
    // determine inside/outsideness. We could try to find a _better_ q which is
    // farther away, but couldn't they all be bad?
    MatrixXV q = BC.row(AJ(0));
    // In a perfect world, it's enough to test a single point.
    double w;

    // winding_number_3 expects colmajor
    const typename DerivedV::Scalar * Vdata;
    Vdata = V.data();
    Matrix<
      typename DerivedV::Scalar,
      DerivedV::RowsAtCompileTime,
      DerivedV::ColsAtCompileTime,
      ColMajor> Vcol;
    if(DerivedV::IsRowMajor)
    {
      // copy to convert to colmajor
      Vcol = V;
      Vdata = Vcol.data();
    }
    winding_number_3(
      Vdata,V.rows(),
      B.data(),B.rows(),
      q.data(),1,&w);
    return fabs(w)>0.5;
  };

  // Reject components which are completely inside other components
  vector<bool> keep(ncc,true);
  size_t nG = 0;
  // This is O( ncc * ncc * m)
  for(size_t id = 0;id<ncc;id++)
  {
    for(size_t oid = 0;oid<ncc;oid++)
    {
      if(id == oid)
      {
        continue;
      }
      const bool inside = is_component_inside_other(V,BC,vG[id],vJ[id],vG[oid]);
#ifdef IGL_OUTER_HULL_DEBUG
      cout<<id<<" is inside "<<oid<<" ? "<<inside<<endl;
#endif
      keep[id] = keep[id] && !inside;
    }
    if(keep[id])
    {
      nG += vJ[id].rows();
    }
  }

  // collect G and J across components
  G.resize(nG,3);
  J.resize(nG,1);
  {
    size_t off = 0;
    for(Index id = 0;id<(Index)ncc;id++)
    {
      if(keep[id])
      {
        assert(vG[id].rows() == vJ[id].rows());
        G.block(off,0,vG[id].rows(),vG[id].cols()) = vG[id];
        J.block(off,0,vJ[id].rows(),vJ[id].cols()) = vJ[id];
        off += vG[id].rows();
      }
    }
  }
}

//.........这里部分代码省略.........
    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 = 
      sIT.rightCols(2).array() /
      (sIT.rightCols(2).colwise()-sIT.col(0)).array();
    MatrixX2S gamma = 
      sIT.rightCols(2).array() /
      (sIT.rightCols(2).colwise()-sIT.col(1)).array();
    vector<Triplet<BCType> > IJV;
    IJV.reserve(m*4*2);
    for(size_t c = 0;c<2;c++)
    {
      for(size_t t = 0;t<(size_t)m;t++)
      {
        IJV.push_back(Triplet<BCType>(0*2*m+c*m+t,  sT(t,0),  lambda(t,c)));
        IJV.push_back(Triplet<BCType>(0*2*m+c*m+t,sT(t,c+2),1-lambda(t,c)));
        IJV.push_back(Triplet<BCType>(1*2*m+c*m+t,  sT(t,1),   gamma(t,c)));
        IJV.push_back(Triplet<BCType>(1*2*m+c*m+t,sT(t,c+2),1- gamma(t,c)));
      }
    }
    BC.resize(m*4,V.rows());
    BC.reserve(m*4*2);
    BC.setFromTriplets(IJV.begin(),IJV.end());
    G.resize(2*m,3);
    G.block(0,0,m,1) = VectorXI::LinSpaced(m,0+0*m,(m-1)+0*m);
    G.block(0,1,m,1) = VectorXI::LinSpaced(m,0+1*m,(m-1)+1*m);
    G.block(0,2,m,1) = VectorXI::LinSpaced(m,0+3*m,(m-1)+3*m);
    G.block(m,0,m,1) = VectorXI::LinSpaced(m,0+0*m,(m-1)+0*m);
    G.block(m,1,m,1) = VectorXI::LinSpaced(m,0+3*m,(m-1)+3*m);
    G.block(m,2,m,1) = VectorXI::LinSpaced(m,0+2*m,(m-1)+2*m);
  };

  MatrixX3I G13,G31,G22;
  SparseMatrix<BCType> BC13,BC31,BC22;
  one_below(T13,IT13,G13,BC13);
  one_below(T31,-IT31,G31,BC31);
  two_below(T22,IT22,G22,BC22);

  BC = cat(1,cat(1,BC13,BC31),BC22);
  U = BC*V;
  G.resize(G13.rows()+G31.rows()+G22.rows(),3);
  G<<G13,(G31.array()+BC13.rows()),(G22.array()+BC13.rows()+BC31.rows());
  MatrixX3S N;
  per_face_normals(U,G,N);
  Matrix<Scalar,1,3> planeN(plane(0),plane(1),plane(2));
  VectorXb flip = (N.array().rowwise() * planeN.array()).rowwise().sum()<0;
  for(size_t g = 0;g<(size_t)G.rows();g++)
  {
    if(flip(g))
    {
      G.row(g) = G.row(g).reverse().eval();
    }
  }

  J.resize(G.rows());
  J<<J13,J31,J22,J22;
}

IGL_INLINE bool igl::biharmonic_coordinates(
  const Eigen::PlainObjectBase<DerivedV> & V,
  const Eigen::PlainObjectBase<DerivedT> & T,
  const std::vector<std::vector<SType> > & S,
  const int k,
  Eigen::PlainObjectBase<DerivedW> & W)
{
  using namespace Eigen;
  using namespace std;
  // This is not the most efficient way to build A, but follows "Linear
  // Subspace Design for Real-Time Shape Deformation" [Wang et al. 2015].
  SparseMatrix<double> A;
  {
    SparseMatrix<double> N,Z,L,K,M;
    normal_derivative(V,T,N);
    Array<bool,Dynamic,1> I;
    Array<bool,Dynamic,Dynamic> C;
    on_boundary(T,I,C);
    {
      std::vector<Triplet<double> >ZIJV;
      for(int t =0;t<T.rows();t++)
      {
        for(int f =0;f<T.cols();f++)
        {
          if(C(t,f))
          {
            const int i = t+f*T.rows();
            for(int c = 1;c<T.cols();c++)
            {
              ZIJV.emplace_back(T(t,(f+c)%T.cols()),i,1);
            }
          }
        }
      }
      Z.resize(V.rows(),N.rows());
      Z.setFromTriplets(ZIJV.begin(),ZIJV.end());
      N = (Z*N).eval();
    }
    cotmatrix(V,T,L);
    K = N+L;
    massmatrix(V,T,MASSMATRIX_TYPE_DEFAULT,M);
    // normalize
    M /= ((VectorXd)M.diagonal()).array().abs().maxCoeff();
    DiagonalMatrix<double,Dynamic> Minv =
      ((VectorXd)M.diagonal().array().inverse()).asDiagonal();
    switch(k)
    {
      default:
        assert(false && "unsupported");
      case 2:
        // For C1 smoothness in 2D, one should use bi-harmonic
        A = K.transpose() * (Minv * K);
        break;
      case 3:
        // For C1 smoothness in 3D, one should use tri-harmonic
        A = K.transpose() * (Minv * (-L * (Minv * K)));
        break;
    }
  }
  // Vertices in point handles
  const size_t mp =
    count_if(S.begin(),S.end(),[](const vector<int> & h){return h.size()==1;});
  // number of region handles
  const size_t r = S.size()-mp;
  // Vertices in region handles
  size_t mr = 0;
  for(const auto & h : S)
  {
    if(h.size() > 1)
    {
      mr += h.size();
    }
  }
  const size_t dim = T.cols()-1;
  // Might as well be dense... I think...
  MatrixXd J = MatrixXd::Zero(mp+mr,mp+r*(dim+1));
  VectorXi b(mp+mr);
  MatrixXd H(mp+r*(dim+1),dim);
  {
    int v = 0;
    int c = 0;
    for(int h = 0;h<S.size();h++)
    {
      if(S[h].size()==1)
      {
        H.row(c) = V.block(S[h][0],0,1,dim);
        J(v,c++) = 1;
        b(v) = S[h][0];
        v++;
      }else
      {
        assert(S[h].size() >= dim+1);
        for(int p = 0;p<S[h].size();p++)
        {
          for(int d = 0;d<dim;d++)
          {
            J(v,c+d) = V(S[h][p],d);
          }
          J(v,c+dim) = 1;
          b(v) = S[h][p];
//.........这里部分代码省略.........