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

IGL_INLINE igl::SolverStatus igl::active_set(
  const Eigen::SparseMatrix<AT>& A,
  const Eigen::PlainObjectBase<DerivedB> & B,
  const Eigen::PlainObjectBase<Derivedknown> & known,
  const Eigen::PlainObjectBase<DerivedY> & Y,
  const Eigen::SparseMatrix<AeqT>& Aeq,
  const Eigen::PlainObjectBase<DerivedBeq> & Beq,
  const Eigen::SparseMatrix<AieqT>& Aieq,
  const Eigen::PlainObjectBase<DerivedBieq> & Bieq,
  const Eigen::PlainObjectBase<Derivedlx> & p_lx,
  const Eigen::PlainObjectBase<Derivedux> & p_ux,
  const igl::active_set_params & params,
  Eigen::PlainObjectBase<DerivedZ> & Z
  )
{
//#define ACTIVE_SET_CPP_DEBUG
#if defined(ACTIVE_SET_CPP_DEBUG) && !defined(_MSC_VER)
#  warning "ACTIVE_SET_CPP_DEBUG"
#endif
  using namespace Eigen;
  using namespace std;
  SolverStatus ret = SOLVER_STATUS_ERROR;
  const int n = A.rows();
  assert(n == A.cols() && "A must be square");
  // Discard const qualifiers
  //if(B.size() == 0)
  //{
  //  B = Eigen::PlainObjectBase<DerivedB>::Zero(n,1);
  //}
  assert(n == B.rows() && "B.rows() must match A.rows()");
  assert(B.cols() == 1 && "B must be a column vector");
  assert(Y.cols() == 1 && "Y must be a column vector");
  assert((Aeq.size() == 0 && Beq.size() == 0) || Aeq.cols() == n);
  assert((Aeq.size() == 0 && Beq.size() == 0) || Aeq.rows() == Beq.rows());
  assert((Aeq.size() == 0 && Beq.size() == 0) || Beq.cols() == 1);
  assert((Aieq.size() == 0 && Bieq.size() == 0) || Aieq.cols() == n);
  assert((Aieq.size() == 0 && Bieq.size() == 0) || Aieq.rows() == Bieq.rows());
  assert((Aieq.size() == 0 && Bieq.size() == 0) || Bieq.cols() == 1);
  Eigen::Matrix<typename Derivedlx::Scalar,Eigen::Dynamic,1> lx;
  Eigen::Matrix<typename Derivedux::Scalar,Eigen::Dynamic,1> ux;
  if(p_lx.size() == 0)
  {
    lx = Eigen::PlainObjectBase<Derivedlx>::Constant(
      n,1,-numeric_limits<typename Derivedlx::Scalar>::max());
  }else
  {
    lx = p_lx;
  }
  if(ux.size() == 0)
  {
    ux = Eigen::PlainObjectBase<Derivedux>::Constant(
      n,1,numeric_limits<typename Derivedux::Scalar>::max());
  }else
  {
    ux = p_ux;
  }
  assert(lx.rows() == n && "lx must have n rows");
  assert(ux.rows() == n && "ux must have n rows");
  assert(ux.cols() == 1 && "lx must be a column vector");
  assert(lx.cols() == 1 && "ux must be a column vector");
  assert((ux.array()-lx.array()).minCoeff() > 0 && "ux(i) must be > lx(i)");
  if(Z.size() != 0)
  {
    // Initial guess should have correct size
    assert(Z.rows() == n && "Z must have n rows");
    assert(Z.cols() == 1 && "Z must be a column vector");
  }
  assert(known.cols() == 1 && "known must be a column vector");
  // Number of knowns
  const int nk = known.size();

  // Initialize active sets
  typedef int BOOL;
#define TRUE 1
#define FALSE 0
  Matrix<BOOL,Dynamic,1> as_lx = Matrix<BOOL,Dynamic,1>::Constant(n,1,FALSE);
  Matrix<BOOL,Dynamic,1> as_ux = Matrix<BOOL,Dynamic,1>::Constant(n,1,FALSE);
  Matrix<BOOL,Dynamic,1> as_ieq = Matrix<BOOL,Dynamic,1>::Constant(Aieq.rows(),1,FALSE);

  // Keep track of previous Z for comparison
  PlainObjectBase<DerivedZ> old_Z;
  old_Z = PlainObjectBase<DerivedZ>::Constant(
      n,1,numeric_limits<typename DerivedZ::Scalar>::max());

  int iter = 0;
  while(true)
  {
#ifdef ACTIVE_SET_CPP_DEBUG
    cout<<"Iteration: "<<iter<<":"<<endl;
    cout<<"  pre"<<endl;
#endif
    // FIND BREACHES OF CONSTRAINTS
    int new_as_lx = 0;
    int new_as_ux = 0;
    int new_as_ieq = 0;
    if(Z.size() > 0)
    {
      for(int z = 0;z < n;z++)
      {
        if(Z(z) < lx(z))
//.........这里部分代码省略.........

IGL_INLINE void igl::orientable_patches(
  const Eigen::PlainObjectBase<DerivedF> & F,
  Eigen::PlainObjectBase<DerivedC> & C,
  Eigen::SparseMatrix<AScalar> & A)
{
  using namespace Eigen;
  using namespace std;

  // simplex size
  assert(F.cols() == 3);

  // List of all "half"-edges: 3*#F by 2
  Matrix<typename DerivedF::Scalar, Dynamic, 2> allE,sortallE,uE;
  allE.resize(F.rows()*3,2);
  Matrix<int,Dynamic,2> IX;
  VectorXi IA,IC;
  allE.block(0*F.rows(),0,F.rows(),1) = F.col(1);
  allE.block(0*F.rows(),1,F.rows(),1) = F.col(2);
  allE.block(1*F.rows(),0,F.rows(),1) = F.col(2);
  allE.block(1*F.rows(),1,F.rows(),1) = F.col(0);
  allE.block(2*F.rows(),0,F.rows(),1) = F.col(0);
  allE.block(2*F.rows(),1,F.rows(),1) = F.col(1);
  // Sort each row
  sort(allE,2,true,sortallE,IX);
  //IC(i) tells us where to find sortallE(i,:) in uE: 
  // so that sortallE(i,:) = uE(IC(i),:)
  unique_rows(sortallE,uE,IA,IC);
  // uE2FT(e,f) = 1 means face f is adjacent to unique edge e
  vector<Triplet<AScalar> > uE2FTijv(IC.rows());
  for(int e = 0;e<IC.rows();e++)
  {
    uE2FTijv[e] = Triplet<AScalar>(e%F.rows(),IC(e),1);
  }
  SparseMatrix<AScalar> uE2FT(F.rows(),uE.rows());
  uE2FT.setFromTriplets(uE2FTijv.begin(),uE2FTijv.end());
  // kill non-manifold edges
  for(int j=0; j<(int)uE2FT.outerSize();j++)
  {
    int degree = 0;
    for(typename SparseMatrix<AScalar>::InnerIterator it (uE2FT,j); it; ++it)
    {
      degree++;
    }
    // Iterate over inside
    if(degree > 2)
    {
      for(typename SparseMatrix<AScalar>::InnerIterator it (uE2FT,j); it; ++it)
      {
        uE2FT.coeffRef(it.row(),it.col()) = 0;
      }
    }
  }
  // Face-face Adjacency matrix
  SparseMatrix<AScalar> uE2F;
  uE2F = uE2FT.transpose().eval();
  A = uE2FT*uE2F;
  // All ones
  for(int j=0; j<A.outerSize();j++)
  {
    // Iterate over inside
    for(typename SparseMatrix<AScalar>::InnerIterator it (A,j); it; ++it)
    {
      if(it.value() > 1)
      {
        A.coeffRef(it.row(),it.col()) = 1;
      }
    }
  }
  //% Connected components are patches
  //%C = components(A); % alternative to graphconncomp from matlab_bgl
  //[~,C] = graphconncomp(A);
  // graph connected components using boost
  components(A,C);

}

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 std;
  using namespace Eigen;

  Eigen::Matrix<
    typename DerivedF::Scalar,Eigen::Dynamic,Eigen::Dynamic>
    FF,FFi;
  triangle_triangle_adjacency(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;

  // 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;
  }

}

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];
//.........这里部分代码省略.........

IGL_INLINE bool igl::bbw::bbw(
  const Eigen::PlainObjectBase<DerivedV> & V,
  const Eigen::PlainObjectBase<DerivedEle> & Ele,
  const Eigen::PlainObjectBase<Derivedb> & b,
  const Eigen::PlainObjectBase<Derivedbc> & bc,
  igl::bbw::BBWData & data,
  Eigen::PlainObjectBase<DerivedW> & W
  )
{
  using namespace std;
  using namespace Eigen;

  // number of domain vertices
  int n = V.rows();
  // number of handles
  int m = bc.cols();
  // Build biharmonic operator
  SparseMatrix<typename DerivedW::Scalar> L;
  cotmatrix(V,Ele,L);
  SparseMatrix<typename DerivedW::Scalar> M;
  SparseMatrix<typename DerivedW::Scalar> Mi;
  massmatrix(V,Ele,MASSMATRIX_TYPE_DEFAULT,M);
  invert_diag(M,Mi);
  SparseMatrix<typename DerivedW::Scalar> Q = L.transpose() * Mi * L;
  assert(!data.partition_unity && "partition_unity not implemented yet");

  W.derived().resize(n,m);
  {
    // No linear terms
    VectorXd c = VectorXd::Zero(n);
    // No linear constraints
    SparseMatrix<typename DerivedW::Scalar> A(0,n),Aeq(0,n),Aieq(0,n);
    VectorXd uc(0,1),Beq(0,1),Bieq(0,1),lc(0,1);
    // Upper and lower box constraints (Constant bounds)
    VectorXd ux = VectorXd::Ones(n);
    VectorXd lx = VectorXd::Zero(n);
    active_set_params eff_params = data.active_set_params;
    switch(data.qp_solver)
    {
      case QP_SOLVER_IGL_ACTIVE_SET:
      {
        if(data.verbosity >= 1)
        {
          cout<<"BBW: max_iter: "<<data.active_set_params.max_iter<<endl;
          cout<<"BBW: eff_max_iter: "<<eff_params.max_iter<<endl;
        }
        if(data.verbosity >= 1)
        {
          cout<<"BBW: Computing initial weights for "<<m<<" handle"<<
            (m!=1?"s":"")<<"."<<endl;
        }
        min_quad_with_fixed_data<typename DerivedW::Scalar > mqwf;
        min_quad_with_fixed_precompute(Q,b,Aeq,true,mqwf);
        min_quad_with_fixed_solve(mqwf,c,bc,Beq,W);
        // decrement
        eff_params.max_iter--;
        bool error = false;
        // Loop over handles
#pragma omp parallel for
        for(int i = 0;i<m;i++)
        {
          // Quicker exit for openmp
          if(error)
          {
            continue;
          }
          if(data.verbosity >= 1)
          {
#pragma omp critical
            cout<<"BBW: Computing weight for handle "<<i+1<<" out of "<<m<<
              "."<<endl;
          }
          VectorXd bci = bc.col(i);
          VectorXd Wi;
          // use initial guess
          Wi = W.col(i);
          SolverStatus ret = active_set(
              Q,c,b,bci,Aeq,Beq,Aieq,Bieq,lx,ux,eff_params,Wi);
          switch(ret)
          {
            case SOLVER_STATUS_CONVERGED:
              break;
            case SOLVER_STATUS_MAX_ITER:
              cerr<<"active_set: max iter without convergence."<<endl;
              break;
            case SOLVER_STATUS_ERROR:
            default:
              cerr<<"active_set error."<<endl;
              error = true;
          }
          W.col(i) = Wi;
        }
        if(error)
        {
          return false;
        }
        break;
      }
      case QP_SOLVER_MOSEK:
      {
//.........这里部分代码省略.........

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

    //igl::sort(di[ui],true,di[ui],IM);
    // Sort, but break ties using index to ensure that duplicates always show
    // up in same order.
    MatrixXd s_di_I;
    igl::sortrows(di_I,true,s_di_I,IM);
    di[ui].resize(uE2E[ui].size());
    for(size_t i = 0;i<di[ui].size();i++)
    {
      di[ui][i] = s_di_I(i,0);
    }

    // copy old list
    vector<typename DerivedF::Index> temp = uE2E[ui];
    for(size_t fei = 0;fei<uE2E[ui].size();fei++)
    {
      uE2E[ui][fei] = temp[IM(fei)];
      const auto & fe = uE2E[ui][fei];
      diIM(fe) = fei;
      dicons(fe) = cons[IM(fei)];
    }

  }

  vector<vector<vector<Index > > > TT,_1;
  triangle_triangle_adjacency(E,EMAP,uE2E,false,TT,_1);
  VectorXI counts;
#ifdef IGL_OUTER_HULL_DEBUG
  cout<<"facet components..."<<endl;
#endif
  facet_components(TT,C,counts);
  assert(C.maxCoeff()+1 == counts.rows());
  const size_t ncc = counts.rows();
  G.resize(0,F.cols());
  J.resize(0,1);
  flip.setConstant(m,1,false);

#ifdef IGL_OUTER_HULL_DEBUG
  cout<<"reindex..."<<endl;
#endif
  // H contains list of faces on outer hull;
  vector<bool> FH(m,false);
  vector<bool> EH(3*m,false);
  vector<MatrixXG> vG(ncc);
  vector<MatrixXJ> vJ(ncc);
  vector<MatrixXJ> vIM(ncc);
  for(size_t id = 0;id<ncc;id++)
  {
    vIM[id].resize(counts[id],1);
  }
  // current index into each IM
  vector<size_t> g(ncc,0);
  // place order of each face in its respective component
  for(Index f = 0;f<m;f++)
  {
    vIM[C(f)](g[C(f)]++) = f;
  }

#ifdef IGL_OUTER_HULL_DEBUG
  cout<<"barycenters..."<<endl;
#endif
  // assumes that "resolve" has handled any coplanar cases correctly and nearly
  // coplanar cases can be sorted based on barycenter.
  MatrixXV BC;
  barycenter(V,F,BC);

IGL_INLINE void igl::bounding_box(
  const Eigen::PlainObjectBase<DerivedV>& V,
  Eigen::PlainObjectBase<DerivedBV>& BV,
  Eigen::PlainObjectBase<DerivedBF>& BF)
{
  using namespace std;

  const int dim = V.cols();
  const auto & minV = V.colwise().minCoeff();
  const auto & maxV = V.colwise().maxCoeff();
  // 2^n vertices
  BV.resize((1<<dim),dim);

  // Recursive lambda to generate all 2^n combinations
  const std::function<void(const int,const int,int*,int)> combos =
  [&BV,&minV,&maxV,&combos](
    const int dim,
    const int i,
    int * X,
    const int pre_index)
  {
    for(X[i] = 0;X[i]<2;X[i]++)
    {
      int index = pre_index*2+X[i];
      if((i+1)<dim)
      {
        combos(dim,i+1,X,index);
      }else
      {
        for(int d = 0;d<dim;d++)
        {
          BV(index,d) = (X[d]?minV[d]:maxV[d]);
        }
      }
    }
  };

  Eigen::VectorXi X(dim);
  combos(dim,0,X.data(),0);
  switch(dim)
  {
    case 2:
      BF.resize(4,2);
      BF<<
        3,1,
        1,0,
        0,2,
        2,3;
      break;
    case 3:
      BF.resize(12,3);
      BF<<
        2,0,6,
        0,4,6,
        5,4,0,
        5,0,1,
        6,4,5,
        5,7,6,
        3,0,2,
        1,0,3,
        3,2,6,
        6,7,3,
        5,1,3,
        3,7,5;
      break;
    default:
      assert(false && "Unsupported dimension.");
      break;
  }
}

IGL_INLINE bool igl::writeOBJ(
  const std::string str,
  const Eigen::PlainObjectBase<DerivedV>& V,
  const Eigen::PlainObjectBase<DerivedF>& F,
  const Eigen::PlainObjectBase<DerivedV>& CN,
  const Eigen::PlainObjectBase<DerivedF>& FN,
  const Eigen::PlainObjectBase<DerivedT>& TC,
  const Eigen::PlainObjectBase<DerivedF>& FTC)
{
  FILE * obj_file = fopen(str.c_str(),"w");
  if(NULL==obj_file)
  {
    printf("IOError: %s could not be opened for writing...",str.c_str());
    return false;
  }
  // Loop over V
  for(int i = 0;i<(int)V.rows();i++)
  {
    fprintf(obj_file,"v %0.15g %0.15g %0.15g\n",
      V(i,0),
      V(i,1),
      V(i,2)
      );
  }
  bool write_N = CN.rows() >0;

  if(write_N)
  {
    for(int i = 0;i<(int)CN.rows();i++)
    {
      fprintf(obj_file,"v %0.15g %0.15g %0.15g\n",
              CN(i,0),
              CN(i,1),
              CN(i,2)
              );
    }
    fprintf(obj_file,"\n");
  }

  bool write_texture_coords = TC.rows() >0;

  if(write_texture_coords)
  {
    for(int i = 0;i<(int)TC.rows();i++)
    {
      fprintf(obj_file, "vt %0.15g %0.15g\n",TC(i,0),TC(i,1));
    }
    fprintf(obj_file,"\n");
  }

  // loop over F
  for(int i = 0;i<(int)F.rows();++i)
  {
    fprintf(obj_file,"f");
    for(int j = 0; j<(int)F.cols();++j)
    {
      // OBJ is 1-indexed
      fprintf(obj_file," %u",F(i,j)+1);

      if(write_texture_coords)
        fprintf(obj_file,"/%u",FTC(i,j)+1);
      if(write_N)
      {
        if (write_texture_coords)
          fprintf(obj_file,"/%u",FN(i,j)+1);
        else
          fprintf(obj_file,"//%u",FN(i,j)+1);
      }
    }
    fprintf(obj_file,"\n");
  }
  fclose(obj_file);
  return true;
}

IGL_INLINE
void igl::copyleft::cgal::order_facets_around_edge(
  const Eigen::PlainObjectBase<DerivedV>& V,
  const Eigen::PlainObjectBase<DerivedF>& F,
  size_t s,
  size_t d, 
  const std::vector<int>& adj_faces,
  const Eigen::PlainObjectBase<DerivedV>& pivot_point,
  Eigen::PlainObjectBase<DerivedI>& order)
{
  assert(V.cols() == 3);
  assert(F.cols() == 3);
  assert(pivot_point.cols() == 3);
  auto signed_index_to_index = [&](int signed_idx)
  {
      return abs(signed_idx) -1;
  };
  auto get_opposite_vertex_index = [&](size_t fid) -> typename DerivedF::Scalar
  {
      typedef typename DerivedF::Scalar Index;
      if (F(fid, 0) != (Index)s && F(fid, 0) != (Index)d) return F(fid, 0);
      if (F(fid, 1) != (Index)s && F(fid, 1) != (Index)d) return F(fid, 1);
      if (F(fid, 2) != (Index)s && F(fid, 2) != (Index)d) return F(fid, 2);
      assert(false);
      // avoid warning
      return -1;
  };

  {
    // Check if s, d and pivot are collinear.
    typedef CGAL::Exact_predicates_exact_constructions_kernel K;
    K::Point_3 ps(V(s,0), V(s,1), V(s,2));
    K::Point_3 pd(V(d,0), V(d,1), V(d,2));
    K::Point_3 pp(pivot_point(0,0), pivot_point(0,1), pivot_point(0,2));
    if (CGAL::collinear(ps, pd, pp)) {
        throw std::runtime_error(
                "Pivot point is collinear with the outer edge!");
    }
  }

  const size_t N = adj_faces.size();
  const size_t num_faces = N + 1; // N adj faces + 1 pivot face

  // Because face indices are used for tie breaking, the original face indices
  // in the new faces array must be ascending.
  auto comp = [&](int i, int j) 
  {
    return signed_index_to_index(adj_faces[i]) <
      signed_index_to_index(adj_faces[j]);
  };
  std::vector<size_t> adj_order(N);
  for (size_t i=0; i<N; i++) adj_order[i] = i;
  std::sort(adj_order.begin(), adj_order.end(), comp);

  DerivedV vertices(num_faces + 2, 3);
  for (size_t i=0; i<N; i++) 
  {
    const size_t fid = signed_index_to_index(adj_faces[adj_order[i]]);
    vertices.row(i) = V.row(get_opposite_vertex_index(fid));
  }
  vertices.row(N  ) = pivot_point;
  vertices.row(N+1) = V.row(s);
  vertices.row(N+2) = V.row(d);

  DerivedF faces(num_faces, 3);
  for (size_t i=0; i<N; i++)
  {
    if (adj_faces[adj_order[i]] < 0) 
    {
      faces(i,0) = N+1; // s
      faces(i,1) = N+2; // d
      faces(i,2) = i  ;
    } else 
    {
      faces(i,0) = N+2; // d
      faces(i,1) = N+1; // s
      faces(i,2) = i  ;
    }
  }
  // Last face is the pivot face.
  faces(N, 0) = N+1;
  faces(N, 1) = N+2;
  faces(N, 2) = N;

  std::vector<int> adj_faces_with_pivot(num_faces);
  for (size_t i=0; i<num_faces; i++)
  {
    if ((size_t)faces(i,0) == N+1 && (size_t)faces(i,1) == N+2)
    {
        adj_faces_with_pivot[i] = int(i+1) * -1;
    } else
    {
        adj_faces_with_pivot[i] = int(i+1);
    }
  }

  DerivedI order_with_pivot;
  order_facets_around_edge(
    vertices, faces, N+1, N+2, adj_faces_with_pivot, order_with_pivot);

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

IGL_INLINE void igl::ismember_rows(
  const Eigen::PlainObjectBase<DerivedA> & A,
  const Eigen::PlainObjectBase<DerivedB> & B,
  Eigen::PlainObjectBase<DerivedIA> & IA,
  Eigen::PlainObjectBase<DerivedLOCB> & LOCB)
{
  using namespace Eigen;
  using namespace std;
  assert(A.cols() == B.cols() && "number of columns must match");
  IA.resize(A.rows(),1);
  IA.setConstant(false);
  LOCB.resize(A.rows(),1);
  LOCB.setConstant(-1);
  // boring base cases
  if(A.size() == 0)
  {
    return;
  }
  if(B.size() == 0)
  {
    return;
  }

  // Get rid of any duplicates
  DerivedA uA;
  DerivedB uB;
  Eigen::Matrix<typename DerivedA::Index,Dynamic,1> uIA,uIuA,uIB,uIuB;
  unique_rows(A,uA,uIA,uIuA);
  unique_rows(B,uB,uIB,uIuB);
  // Sort both
  DerivedA sA;
  DerivedB sB;
  Eigen::Matrix<typename DerivedA::Index,Dynamic,1> sIA,sIB;
  sortrows(uA,true,sA,sIA);
  sortrows(uB,true,sB,sIB);

  Eigen::Matrix<bool,Eigen::Dynamic,1> uF = 
    Eigen::Matrix<bool,Eigen::Dynamic,1>::Zero(sA.size(),1);
  Eigen::Matrix<typename DerivedLOCB::Scalar, Eigen::Dynamic,1> uLOCB =
    Eigen::Matrix<typename DerivedLOCB::Scalar,Eigen::Dynamic,1>::
    Constant(sA.size(),1,-1);
  const auto & row_greater_than = [&sA,&sB](const int a, const int b)
  {
    for(int c = 0;c<sA.cols();c++)
    {
      if(sA(a,c) > sB(b,c)) return true;
      if(sA(a,c) < sB(b,c)) return false;
    }
    return false;
  };
  {
    int bi = 0;
    // loop over sA
    bool past = false;
    for(int a = 0;a<sA.rows();a++)
    {
      while(!past && row_greater_than(a,bi))
      {
        bi++;
        past = bi>=sB.size();
      }
      if(!past && (sA.row(a).array()==sB.row(bi).array()).all() )
      {
        uF(sIA(a)) = true;
        uLOCB(sIA(a)) = uIB(sIB(bi));
      }
    }
  }

  for(int a = 0;a<A.rows();a++)
  {
    IA(a) = uF(uIuA(a));
    LOCB(a) = uLOCB(uIuA(a));
  }
}

IGL_INLINE void igl::ismember(
  const Eigen::MatrixBase<DerivedA> & A,
  const Eigen::MatrixBase<DerivedB> & B,
  Eigen::PlainObjectBase<DerivedIA> & IA,
  Eigen::PlainObjectBase<DerivedLOCB> & LOCB)
{
  using namespace Eigen;
  using namespace std;
  IA.resizeLike(A);
  IA.setConstant(false);
  LOCB.resizeLike(A);
  LOCB.setConstant(-1);
  // boring base cases
  if(A.size() == 0)
  {
    return;
  }
  if(B.size() == 0)
  {
    return;
  }

  // Get rid of any duplicates
  typedef Matrix<typename DerivedA::Scalar,Dynamic,1> VectorA;
  typedef Matrix<typename DerivedB::Scalar,Dynamic,1> VectorB;
  const VectorA vA(Eigen::Map<const VectorA>(DerivedA(A).data(), A.cols()*A.rows(),1));
  const VectorB vB(Eigen::Map<const VectorB>(DerivedB(B).data(), B.cols()*B.rows(),1));
  VectorA uA;
  VectorB uB;
  Eigen::Matrix<typename DerivedA::Index,Dynamic,1> uIA,uIuA,uIB,uIuB;
  unique(vA,uA,uIA,uIuA);
  unique(vB,uB,uIB,uIuB);
  // Sort both
  VectorA sA;
  VectorB sB;
  Eigen::Matrix<typename DerivedA::Index,Dynamic,1> sIA,sIB;
  sort(uA,1,true,sA,sIA);
  sort(uB,1,true,sB,sIB);

  Eigen::Matrix<bool,Eigen::Dynamic,1> uF = 
    Eigen::Matrix<bool,Eigen::Dynamic,1>::Zero(sA.size(),1);
  Eigen::Matrix<typename DerivedLOCB::Scalar, Eigen::Dynamic,1> uLOCB =
    Eigen::Matrix<typename DerivedLOCB::Scalar,Eigen::Dynamic,1>::
    Constant(sA.size(),1,-1);
  {
    int bi = 0;
    // loop over sA
    bool past = false;
    for(int a = 0;a<sA.size();a++)
    {
      while(!past && sA(a)>sB(bi))
      {
        bi++;
        past = bi>=sB.size();
      }
      if(!past && sA(a)==sB(bi))
      {
        uF(sIA(a)) = true;
        uLOCB(sIA(a)) = uIB(sIB(bi));
      }
    }
  }

  Map< Matrix<typename DerivedIA::Scalar,Dynamic,1> > 
    vIA(IA.data(),IA.cols()*IA.rows(),1);
  Map< Matrix<typename DerivedLOCB::Scalar,Dynamic,1> > 
    vLOCB(LOCB.data(),LOCB.cols()*LOCB.rows(),1);
  for(int a = 0;a<A.size();a++)
  {
    vIA(a) = uF(uIuA(a));
    vLOCB(a) = uLOCB(uIuA(a));
  }
}

IGL_INLINE void igl::boolean::mesh_boolean(
  const Eigen::PlainObjectBase<DerivedVA > & VA,
  const Eigen::PlainObjectBase<DerivedFA > & FA,
  const Eigen::PlainObjectBase<DerivedVB > & VB,
  const Eigen::PlainObjectBase<DerivedFB > & FB,
  const MeshBooleanType & type,
  const std::function<void(
    const Eigen::Matrix<typename DerivedVC::Scalar,Eigen::Dynamic,3>&,
    const Eigen::Matrix<typename DerivedFC::Scalar, Eigen::Dynamic,3>&,
          Eigen::Matrix<typename DerivedVC::Scalar,Eigen::Dynamic,3>&,
          Eigen::Matrix<typename DerivedFC::Scalar, Eigen::Dynamic,3>&,
          Eigen::Matrix<typename DerivedJ::Scalar, Eigen::Dynamic,1>&)>
    & resolve_fun,
  Eigen::PlainObjectBase<DerivedVC > & VC,
  Eigen::PlainObjectBase<DerivedFC > & FC,
  Eigen::PlainObjectBase<DerivedJ > & J)
{
  using namespace Eigen;
  using namespace std;
  using namespace igl;
  using namespace igl::cgal;
  MeshBooleanType eff_type = type;
  // Concatenate A and B into a single mesh
  typedef CGAL::Exact_predicates_exact_constructions_kernel Kernel;
  typedef Kernel::FT ExactScalar;
  typedef typename DerivedVC::Scalar Scalar;
  typedef typename DerivedFC::Scalar Index;
  typedef Matrix<Scalar,Dynamic,3> MatrixX3S;
  typedef Matrix<ExactScalar,Dynamic,3> MatrixX3ES;
  typedef Matrix<Index,Dynamic,3> MatrixX3I;
  typedef Matrix<Index,Dynamic,2> MatrixX2I;
  typedef Matrix<Index,Dynamic,1> VectorXI;
  typedef Matrix<typename DerivedJ::Scalar,Dynamic,1> VectorXJ;
#ifdef IGL_MESH_BOOLEAN_DEBUG
  cout<<"mesh boolean..."<<endl;
#endif
  MatrixX3S V(VA.rows()+VB.rows(),3);
  MatrixX3I F(FA.rows()+FB.rows(),3);
  V.block(0,0,VA.rows(),VA.cols()) = VA;
  V.block(VA.rows(),0,VB.rows(),VB.cols()) = VB;
#ifdef IGL_MESH_BOOLEAN_DEBUG
  cout<<"prepare selfintersect input..."<<endl;
#endif
  switch(type)
  {
    // Minus is implemented by flipping B and computing union
    case MESH_BOOLEAN_TYPE_MINUS:
      F.block(0,0,FA.rows(),FA.cols()) = FA.rowwise().reverse();
      F.block(FA.rows(),0,FB.rows(),FB.cols()) = FB.array()+VA.rows();
      //F.block(0,0,FA.rows(),3) = FA;
      //F.block(FA.rows(),0,FB.rows(),3) =
      //  FB.rowwise().reverse().array()+VA.rows();
      eff_type = MESH_BOOLEAN_TYPE_INTERSECT;
      break;
    default:
      F.block(0,0,FA.rows(),FA.cols()) = FA;
      F.block(FA.rows(),0,FB.rows(),FB.cols()) = FB.array()+VA.rows();
      break;
  }

  // Resolve intersections (assumes A and B are solid)
  const auto & libigl_resolve = [](
    const MatrixX3S & V,
    const MatrixX3I & F,
    MatrixX3ES & CV,
    MatrixX3I & CF,
    VectorXJ & J)
  {
    MatrixX3ES SV;
    MatrixX3I SF;
    MatrixX2I SIF;
    VectorXI SIM,UIM;
    igl::cgal::RemeshSelfIntersectionsParam params;
    remesh_self_intersections(V,F,params,SV,SF,SIF,J,SIM);
    for_each(SF.data(),SF.data()+SF.size(),[&SIM](int & a){a=SIM(a);});
    {
      remove_unreferenced(SV,SF,CV,CF,UIM);
    }
  };

#ifdef IGL_MESH_BOOLEAN_DEBUG
  cout<<"resolve..."<<endl;
#endif
  MatrixX3S CV;
  MatrixX3ES EV;
  MatrixX3I CF;
  VectorXJ CJ;
  if(resolve_fun)
  {
    resolve_fun(V,F,CV,CF,CJ);
  }else
  {
    libigl_resolve(V,F,EV,CF,CJ);
    CV.resize(EV.rows(), EV.cols());
    std::transform(EV.data(), EV.data() + EV.rows()*EV.cols(),
            CV.data(), [&](ExactScalar val) {
            return CGAL::to_double(val);
            });
  }

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

IGL_INLINE void igl::orient_outward(
  const Eigen::PlainObjectBase<DerivedV> & V,
  const Eigen::PlainObjectBase<DerivedF> & F,
  const Eigen::PlainObjectBase<DerivedC> & C,
  Eigen::PlainObjectBase<DerivedFF> & FF,
  Eigen::PlainObjectBase<DerivedI> & I)
{
  using namespace Eigen;
  using namespace std;
  assert(C.rows() == F.rows());
  assert(F.cols() == 3);
  assert(V.cols() == 3);

  // number of faces
  const int m = F.rows();
  // number of patches
  const int num_cc = C.maxCoeff()+1;
  I.resize(num_cc);
  if(&FF != &F)
  {
    FF = F;
  }
  PlainObjectBase<DerivedV> N,BC,BCmean;
  Matrix<typename DerivedV::Scalar,Dynamic,1> A;
  VectorXd totA(num_cc), dot(num_cc);
  Matrix<typename DerivedV::Scalar,3,1> Z(1,1,1);
  per_face_normals(V,F,Z.normalized(),N);
  barycenter(V,F,BC);
  doublearea(V,F,A);
  BCmean.setConstant(num_cc,3,0);
  dot.setConstant(num_cc,1,0);
  totA.setConstant(num_cc,1,0);
  // loop over faces
  for(int f = 0;f<m;f++)
  {
    BCmean.row(C(f)) += A(f)*BC.row(f);
    totA(C(f))+=A(f);
  }
  // take area weighted average
  for(int c = 0;c<num_cc;c++)
  {
    BCmean.row(c) /= (typename DerivedV::Scalar) totA(c);
  }
  // subtract bcmean
  for(int f = 0;f<m;f++)
  {
    BC.row(f) -= BCmean.row(C(f));
    dot(C(f)) += A(f)*N.row(f).dot(BC.row(f));
  }
  // take area weighted average
  for(int c = 0;c<num_cc;c++)
  {
    dot(c) /= (typename DerivedV::Scalar) totA(c);
    if(dot(c) < 0)
    {
      I(c) = true;
    }else
    {
      I(c) = false;
    }
  }
  // flip according to I
  for(int f = 0;f<m;f++)
  {
    if(I(C(f)))
    {
      FF.row(f) = FF.row(f).reverse().eval();
    }
  }
}

IGL_INLINE void igl::cotangent(
  const Eigen::PlainObjectBase<DerivedV>& V,
  const Eigen::PlainObjectBase<DerivedF>& F,
  Eigen::PlainObjectBase<DerivedC>& C)
{
  using namespace igl;
  using namespace std;
  using namespace Eigen;
  // simplex size (3: triangles, 4: tetrahedra)
  int simplex_size = F.cols();
  // Number of elements
  int m = F.rows();

  // Law of cosines + law of sines
  switch(simplex_size)
  {
    case 3:
    {
      // Triangles
      //Matrix<typename DerivedC::Scalar,Dynamic,3> l;
      //edge_lengths(V,F,l);
      // edge lengths numbered same as opposite vertices
      Matrix<typename DerivedC::Scalar,Dynamic,3> l;
      igl::edge_lengths(V,F,l);
      // double area
      Matrix<typename DerivedC::Scalar,Dynamic,1> dblA;
      doublearea(l,dblA);
      // cotangents and diagonal entries for element matrices
      // correctly divided by 4 (alec 2010)
      C.resize(m,3);
      for(int i = 0;i<m;i++)
      {
        C(i,0) = (l(i,1)*l(i,1) + l(i,2)*l(i,2) - l(i,0)*l(i,0))/dblA(i)/4.0;
        C(i,1) = (l(i,2)*l(i,2) + l(i,0)*l(i,0) - l(i,1)*l(i,1))/dblA(i)/4.0;
        C(i,2) = (l(i,0)*l(i,0) + l(i,1)*l(i,1) - l(i,2)*l(i,2))/dblA(i)/4.0;
      }
      break;
    }
    case 4:
    {

      // edge lengths numbered same as opposite vertices
      Matrix<typename DerivedC::Scalar,Dynamic,6> l;
      edge_lengths(V,F,l);
      Matrix<typename DerivedC::Scalar,Dynamic,4> s;
      face_areas(l,s);
      Matrix<typename DerivedC::Scalar,Dynamic,6> cos_theta,theta;
      dihedral_angles_intrinsic(l,s,theta,cos_theta);

      // volume
      Matrix<typename DerivedC::Scalar,Dynamic,1> vol;
      volume(l,vol);


      // Law of sines
      // http://mathworld.wolfram.com/Tetrahedron.html
      Matrix<typename DerivedC::Scalar,Dynamic,6> sin_theta(m,6);
      sin_theta.col(0) = vol.array() / ((2./(3.*l.col(0).array())).array() * s.col(1).array() * s.col(2).array());
      sin_theta.col(1) = vol.array() / ((2./(3.*l.col(1).array())).array() * s.col(2).array() * s.col(0).array());
      sin_theta.col(2) = vol.array() / ((2./(3.*l.col(2).array())).array() * s.col(0).array() * s.col(1).array());
      sin_theta.col(3) = vol.array() / ((2./(3.*l.col(3).array())).array() * s.col(3).array() * s.col(0).array());
      sin_theta.col(4) = vol.array() / ((2./(3.*l.col(4).array())).array() * s.col(3).array() * s.col(1).array());
      sin_theta.col(5) = vol.array() / ((2./(3.*l.col(5).array())).array() * s.col(3).array() * s.col(2).array());


      // http://arxiv.org/pdf/1208.0354.pdf Page 18
      C = (1./6.) * l.array() * cos_theta.array() / sin_theta.array();

      break;
    }
    default:
    {
      fprintf(stderr,
          "cotangent.h: Error: Simplex size (%d) not supported\n", simplex_size);
      assert(false);
    }
  }
}

IGL_INLINE void igl::boundary_loop(
    const Eigen::PlainObjectBase<DerivedF> & F, 
    std::vector<std::vector<Index> >& L)
{
  using namespace std;
  using namespace Eigen;
  using namespace igl;

  MatrixXd Vdummy(F.maxCoeff(),1);
  MatrixXi TT,TTi;
  vector<std::vector<int> > VF, VFi;
  triangle_triangle_adjacency(Vdummy,F,TT,TTi);
  vertex_triangle_adjacency(Vdummy,F,VF,VFi);

  vector<bool> unvisited = is_border_vertex(Vdummy,F);
  set<int> unseen;
  for (int i = 0; i < unvisited.size(); ++i)
  {
    if (unvisited[i])
      unseen.insert(unseen.end(),i);
  }

  while (!unseen.empty())
  {
    vector<Index> l;

    // Get first vertex of loop
    int start = *unseen.begin();
    unseen.erase(unseen.begin());
    unvisited[start] = false;
    l.push_back(start);

    bool done = false;
    while (!done)
    {
      // Find next vertex
      bool newBndEdge = false;
      int v = l[l.size()-1];
      int next;
      for (int i = 0; i < (int)VF[v].size() && !newBndEdge; i++)
      {
        int fid = VF[v][i];

        if (TT.row(fid).minCoeff() < 0.) // Face contains boundary edge
        {
          int vLoc;
          if (F(fid,0) == v) vLoc = 0;
          if (F(fid,1) == v) vLoc = 1;
          if (F(fid,2) == v) vLoc = 2;

          int vPrev = F(fid,(vLoc + F.cols()-1) % F.cols());
          int vNext = F(fid,(vLoc + 1) % F.cols());

          bool newBndEdge = false;
          if (unvisited[vPrev] && TT(fid,(vLoc+2) % F.cols()) < 0)
          {
            next = vPrev;
            newBndEdge = true;
          }
          else if (unvisited[vNext] && TT(fid,vLoc) < 0)
          {
            next = vNext;
            newBndEdge = true;
          }
        }
      }

      if (newBndEdge)
      {
        l.push_back(next);
        unseen.erase(next);
        unvisited[next] = false;
      }
      else
        done = true;
    }
    L.push_back(l);
  }
}