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

IGL_INLINE Eigen::PlainObjectBase<DerivedV> igl::dot_row(
  const Eigen::PlainObjectBase<DerivedV>& A,
  const Eigen::PlainObjectBase<DerivedV>& B
  )
{
  assert(A.rows() == B.rows());
  assert(A.cols() == B.cols());

  return (A.array() * B.array()).rowwise().sum();
}

IGL_INLINE void igl::parse_rhs_index(
    const mxArray *prhs[], 
    Eigen::PlainObjectBase<DerivedV> & V)
{
  parse_rhs_double(prhs,V);
  V.array() -= 1;
}

IGL_INLINE void igl::barycentric_coordinates(
  const Eigen::PlainObjectBase<DerivedP> & P,
  const Eigen::PlainObjectBase<DerivedA> & A,
  const Eigen::PlainObjectBase<DerivedB> & B,
  const Eigen::PlainObjectBase<DerivedC> & C,
  const Eigen::PlainObjectBase<DerivedD> & D,
  Eigen::PlainObjectBase<DerivedL> & L)
{
  using namespace Eigen;
  assert(P.cols() == 3 && "query must be in 3d");
  assert(A.cols() == 3 && "corners must be in 3d");
  assert(B.cols() == 3 && "corners must be in 3d");
  assert(C.cols() == 3 && "corners must be in 3d");
  assert(D.cols() == 3 && "corners must be in 3d");
  assert(P.rows() == A.rows() && "Must have same number of queries as corners");
  assert(A.rows() == B.rows() && "Corners must be same size");
  assert(A.rows() == C.rows() && "Corners must be same size");
  assert(A.rows() == D.rows() && "Corners must be same size");
  typedef Matrix<typename DerivedL::Scalar,DerivedL::RowsAtCompileTime,1> 
    VectorXS;
  // Total volume
  VectorXS vol,LA,LB,LC,LD;
  volume(B,D,C,P,LA);
  volume(A,C,D,P,LB);
  volume(A,D,B,P,LC);
  volume(A,B,C,P,LD);
  volume(A,B,C,D,vol);
  L.resize(P.rows(),4);
  L<<LA,LB,LC,LD;
  L.array().colwise() /= vol.array();
}

IGL_INLINE void igl::doublearea(
  const Eigen::PlainObjectBase<DerivedV> & V,
  const Eigen::PlainObjectBase<DerivedF> & F,
  Eigen::PlainObjectBase<DeriveddblA> & dblA)
{
  if (F.cols() == 4) // quads are handled by a specialized function
    return doublearea_quad(V,F,dblA);

  const int dim = V.cols();
  // Only support triangles
  assert(F.cols() == 3);
  const size_t m = F.rows();
  // Compute edge lengths
  Eigen::PlainObjectBase<DerivedV> l;
  // "Lecture Notes on Geometric Robustness" Shewchuck 09, Section 3.1
  // http://www.cs.berkeley.edu/~jrs/meshpapers/robnotes.pdf

  // Projected area helper
  const auto & proj_doublearea =
    [&V,&F](const int x, const int y, const int f)->double
  {
    auto rx = V(F(f,0),x)-V(F(f,2),x);
    auto sx = V(F(f,1),x)-V(F(f,2),x);
    auto ry = V(F(f,0),y)-V(F(f,2),y);
    auto sy = V(F(f,1),y)-V(F(f,2),y);
    return rx*sy - ry*sx;
  };

  switch(dim)
  {
    case 3:
    {
      dblA = Eigen::PlainObjectBase<DeriveddblA>::Zero(m,1);
      for(size_t f = 0;f<m;f++)
      {
        for(int d = 0;d<3;d++)
        {
          double dblAd = proj_doublearea(d,(d+1)%3,f);
          dblA(f) += dblAd*dblAd;
        }
      }
      dblA = dblA.array().sqrt().eval();
      break;
    }
    case 2:
    {
      dblA.resize(m,1);
      for(size_t f = 0;f<m;f++)
      {
        dblA(f) = proj_doublearea(0,1,f);
      }
      break;
    }
    default:
    {
      edge_lengths(V,F,l);
      return doublearea(l,dblA);
    }
  }
}

IGL_INLINE void igl::doublearea( 
  const Eigen::PlainObjectBase<DerivedV> & V, 
  const Eigen::PlainObjectBase<DerivedF> & F, 
  Eigen::PlainObjectBase<DeriveddblA> & dblA)
{
  const int dim = V.cols();
  // Only support triangles
  assert(F.cols() == 3);
  const int m = F.rows();
  // Compute edge lengths
  Eigen::PlainObjectBase<DerivedV> l;
  // "Lecture Notes on Geometric Robustness" Shewchuck 09, Section 3.1
  // http://www.cs.berkeley.edu/~jrs/meshpapers/robnotes.pdf
  switch(dim)
  {
    case 3:
    {
      dblA = Eigen::PlainObjectBase<DeriveddblA>::Zero(m,1);
      for(int d = 0;d<3;d++)
      {
        Eigen::Matrix<typename DerivedV::Scalar, DerivedV::RowsAtCompileTime, 2> Vd(V.rows(),2);
        Vd << V.col(d), V.col((d+1)%3);
        Eigen::PlainObjectBase<DeriveddblA> dblAd;
        doublearea(Vd,F,dblAd);
        dblA.array() += dblAd.array().square();
      }
      dblA = dblA.array().sqrt().eval();
      break;
    }
    case 2:
    {
      dblA.resize(m,1);
      for(int f = 0;f<m;f++)
      {
        auto r = V.row(F(f,0))-V.row(F(f,2));
        auto s = V.row(F(f,1))-V.row(F(f,2));
        dblA(f) = r(0)*s(1) - r(1)*s(0);
      }
      break;
    }
    default:
    {
      edge_lengths(V,F,l);
      return doublearea(l,dblA);
    }
  }
}

IGL_INLINE void igl::dihedral_angles_intrinsic(
  const Eigen::PlainObjectBase<DerivedL>& L,
  const Eigen::PlainObjectBase<DerivedA>& A,
  Eigen::PlainObjectBase<Derivedtheta>& theta,
  Eigen::PlainObjectBase<Derivedcos_theta>& cos_theta)
{
  using namespace Eigen;
  const int m = L.rows();
  assert(m == A.rows());
  // Law of cosines
  // http://math.stackexchange.com/a/49340/35376
  Matrix<typename Derivedtheta::Scalar,Dynamic,6> H_sqr(m,6);
  H_sqr.col(0) = (1./16.) * (4. * L.col(3).array().square() * L.col(0).array().square() - 
                                ((L.col(1).array().square() + L.col(4).array().square()) -
                                 (L.col(2).array().square() + L.col(5).array().square())).square());
  H_sqr.col(1) = (1./16.) * (4. * L.col(4).array().square() * L.col(1).array().square() - 
                                ((L.col(2).array().square() + L.col(5).array().square()) -
                                 (L.col(3).array().square() + L.col(0).array().square())).square());
  H_sqr.col(2) = (1./16.) * (4. * L.col(5).array().square() * L.col(2).array().square() - 
                                ((L.col(3).array().square() + L.col(0).array().square()) -
                                 (L.col(4).array().square() + L.col(1).array().square())).square());
  H_sqr.col(3) = (1./16.) * (4. * L.col(0).array().square() * L.col(3).array().square() - 
                                ((L.col(4).array().square() + L.col(1).array().square()) -
                                 (L.col(5).array().square() + L.col(2).array().square())).square());
  H_sqr.col(4) = (1./16.) * (4. * L.col(1).array().square() * L.col(4).array().square() - 
                                ((L.col(5).array().square() + L.col(2).array().square()) -
                                 (L.col(0).array().square() + L.col(3).array().square())).square());
  H_sqr.col(5) = (1./16.) * (4. * L.col(2).array().square() * L.col(5).array().square() - 
                                ((L.col(0).array().square() + L.col(3).array().square()) -
                                 (L.col(1).array().square() + L.col(4).array().square())).square());
  cos_theta.resize(m,6);
  cos_theta.col(0) = (H_sqr.col(0).array() - 
      A.col(1).array().square() - A.col(2).array().square()).array() / 
      (-2.*A.col(1).array() * A.col(2).array());
  cos_theta.col(1) = (H_sqr.col(1).array() - 
      A.col(2).array().square() - A.col(0).array().square()).array() / 
      (-2.*A.col(2).array() * A.col(0).array());
  cos_theta.col(2) = (H_sqr.col(2).array() - 
      A.col(0).array().square() - A.col(1).array().square()).array() / 
      (-2.*A.col(0).array() * A.col(1).array());
  cos_theta.col(3) = (H_sqr.col(3).array() - 
      A.col(3).array().square() - A.col(0).array().square()).array() / 
      (-2.*A.col(3).array() * A.col(0).array());
  cos_theta.col(4) = (H_sqr.col(4).array() - 
      A.col(3).array().square() - A.col(1).array().square()).array() / 
      (-2.*A.col(3).array() * A.col(1).array());
  cos_theta.col(5) = (H_sqr.col(5).array() - 
      A.col(3).array().square() - A.col(2).array().square()).array() / 
      (-2.*A.col(3).array() * A.col(2).array());

  theta = cos_theta.array().acos();

  cos_theta.resize(m,6);

}

IGL_INLINE bool igl::readBF(
    const std::string & filename,
    Eigen::PlainObjectBase<DerivedWI> & WI,
    Eigen::PlainObjectBase<DerivedbfP> & bfP,
    Eigen::PlainObjectBase<DerivedO> & offsets,
    Eigen::PlainObjectBase<DerivedC> & C,
    Eigen::PlainObjectBase<DerivedBE> & BE,
    Eigen::PlainObjectBase<DerivedP> & P)
{
    using namespace Eigen;
    using namespace std;
    if(!readBF(filename,WI,bfP,offsets))
    {
        return false;
    }

    C.resize(WI.rows(),3);
    vector<bool> computed(C.rows(),false);
    // better not be cycles in bfP
    std::function<Eigen::RowVector3d(const int)> locate_tip;
    locate_tip =
        [&offsets,&computed,&bfP,&locate_tip,&C](const int w)->Eigen::RowVector3d
    {
        if(w<0) return Eigen::RowVector3d(0,0,0);
        if(computed[w]) return C.row(w);
        computed[w] = true;
        return C.row(w) = locate_tip(bfP(w)) + offsets.row(w);
    };
    int num_roots = (bfP.array() == -1).count();
    BE.resize(WI.rows()-num_roots,2);
    P.resize(BE.rows());
    for(int c = 0; c<C.rows(); c++)
    {
        locate_tip(c);
        assert(c>=0);
        // weight associated with this bone
        const int wi = WI(c);
        if(wi >= 0)
        {
            // index into C
            const int p = bfP(c);
            assert(p >= 0 && "No weights for roots allowed");
            // index into BE
            const int pwi = WI(p);
            P(wi) = pwi;
            BE(wi,0) = p;
            BE(wi,1) = c;
        }
    }
    return true;
}

IGL_INLINE bool igl::writeOBJ(
  const std::string str,
  const Eigen::PlainObjectBase<DerivedV>& V,
  const Eigen::PlainObjectBase<DerivedF>& F)
{
  using namespace std;
  using namespace Eigen;
  assert(V.cols() == 3 && "V should have 3 columns");
  ofstream s(str);
  if(!s.is_open())
  {
    fprintf(stderr,"IOError: writeOBJ() could not open %s\n",str.c_str());
    return false;
  }
  s<<
    V.format(IOFormat(FullPrecision,DontAlignCols," ","\n","v ","","","\n"))<<
    (F.array()+1).format(IOFormat(FullPrecision,DontAlignCols," ","\n","f ","","","\n"));
  return true;
}

void igl::copyleft::offset_surface(
  const Eigen::MatrixBase<DerivedV> & V,
  const Eigen::MatrixBase<DerivedF> & F,
  const isolevelType isolevel,
  const typename Derivedside::Scalar s,
  const SignedDistanceType & signed_distance_type,
  Eigen::PlainObjectBase<DerivedSV> & SV,
  Eigen::PlainObjectBase<DerivedSF> & SF,
  Eigen::PlainObjectBase<DerivedGV> &   GV,
  Eigen::PlainObjectBase<Derivedside> & side,
  Eigen::PlainObjectBase<DerivedS> & S)
{
  typedef typename DerivedV::Scalar Scalar;
  typedef typename DerivedF::Scalar Index;
  {
    Eigen::AlignedBox<Scalar,3> box;
    typedef Eigen::Matrix<Scalar,1,3> RowVector3S;
    assert(V.cols() == 3 && "V must contain positions in 3D");
    RowVector3S min_ext = V.colwise().minCoeff().array() - isolevel;
    RowVector3S max_ext = V.colwise().maxCoeff().array() + isolevel;
    box.extend(min_ext.transpose());
    box.extend(max_ext.transpose());
    igl::voxel_grid(box,s,1,GV,side);
  }

  const Scalar h = 
    (GV.col(0).maxCoeff()-GV.col(0).minCoeff())/((Scalar)(side(0)-1));
  const Scalar lower_bound = isolevel-sqrt(3.0)*h;
  const Scalar upper_bound = isolevel+sqrt(3.0)*h;
  {
    Eigen::Matrix<Index,Eigen::Dynamic,1> I;
    Eigen::Matrix<typename DerivedV::Scalar,Eigen::Dynamic,3> C,N;
    igl::signed_distance(
      GV,V,F,signed_distance_type,lower_bound,upper_bound,S,I,C,N);
  }
  igl::flood_fill(side,S);
  
  DerivedS SS = S.array()-isolevel;
  igl::copyleft::marching_cubes(SS,GV,side(0),side(1),side(2),SV,SF);
}

//.........这里部分代码省略.........
        OffendingMap offendingA,offendingB;
        //EdgeMap edge2facesA,edge2facesB;

        std::list<int> lIF;
        const auto cb = [&](const Box &a, const Box &b) -> void
        {
          using namespace std;
          // index in F and T
          int fa = a.handle()-TA.begin();
          int fb = b.handle()-TB.begin();
          const Triangle_3 & A = *a.handle();
          const Triangle_3 & B = *b.handle();
          if(CGAL::do_intersect(A,B))
          {
            // There was an intersection
            lIF.push_back(fa);
            lIF.push_back(fb);
            if(params.first_only)
            {
              throw IGL_FIRST_HIT_EXCEPTION;
            }
            if(!params.detect_only)
            {
              CGAL::Object result = CGAL::intersection(A,B);

              push_result(FA,fa,fb,result,offendingA);
              push_result(FB,fb,fa,result,offendingB);
            }
          }
        };
        try{
          CGAL::box_intersection_d(
            A_boxes.begin(), A_boxes.end(),
            B_boxes.begin(), B_boxes.end(),
            cb);
        }catch(int e)
        {
          // Rethrow if not FIRST_HIT_EXCEPTION
          if(e != IGL_FIRST_HIT_EXCEPTION)
          {
            throw e;
          }
          // Otherwise just fall through
        }

        // Convert lIF to Eigen matrix
        assert(lIF.size()%2 == 0);
        IF.resize(lIF.size()/2,2);
        {
          int i=0;
          for(
            list<int>::const_iterator ifit = lIF.begin();
            ifit!=lIF.end();
            )
          {
            IF(i,0) = (*ifit);
            ifit++; 
            IF(i,1) = (*ifit);
            ifit++;
            i++;
          }
        }
        if(!params.detect_only)
        {
          // Obsolete, now remesh_intersections expects a single mesh
          // remesh_intersections(VA,FA,TA,offendingA,VVA,FFA,JA,IMA);
          // remesh_intersections(VB,FB,TB,offendingB,VVB,FFB,JB,IMB);
          // Combine mesh and offending maps
          DerivedVA VAB(VA.rows()+VB.rows(),3);
          VAB<<VA,VB;
          DerivedFA FAB(FA.rows()+FB.rows(),3);
          FAB<<FA,(FB.array()+VA.rows());
          Triangles TAB;
          TAB.reserve(TA.size()+TB.size());
          TAB.insert(TAB.end(),TA.begin(),TA.end());
          TAB.insert(TAB.end(),TB.begin(),TB.end());
          OffendingMap offending;
          //offending.reserve(offendingA.size() + offendingB.size());
          for (const auto itr : offendingA)
          {
            // Remap offenders in FB to FAB
            auto offenders = itr.second;
            for(auto & offender : offenders)
            {
              offender.first += FA.rows();
            }
            offending[itr.first] = offenders;
          }
          for (const auto itr : offendingB)
          {
            // Store offenders for FB according to place in FAB
            offending[FA.rows() + itr.first] = itr.second;
          }

          remesh_intersections(
            VAB,FAB,TAB,offending,params.stitch_all,VVAB,FFAB,JAB,IMAB);
        }

        return IF.rows() > 0;
      }

IGL_INLINE bool igl::copyleft::boolean::mesh_boolean(
    const Eigen::PlainObjectBase<DerivedVA> & VA,
    const Eigen::PlainObjectBase<DerivedFA> & FA,
    const Eigen::PlainObjectBase<DerivedVB> & VB,
    const Eigen::PlainObjectBase<DerivedFB> & FB,
    const WindingNumberOp& wind_num_op,
    const KeepFunc& keep,
    const ResolveFunc& resolve_fun,
    Eigen::PlainObjectBase<DerivedVC > & VC,
    Eigen::PlainObjectBase<DerivedFC > & FC,
    Eigen::PlainObjectBase<DerivedJ > & J) 
{

#ifdef MESH_BOOLEAN_TIMING
  const auto & tictoc = []() -> double
  {
    static double t_start = igl::get_seconds();
    double diff = igl::get_seconds()-t_start;
    t_start += diff;
    return diff;
  };
  const auto log_time = [&](const std::string& label) -> void {
    std::cout << "mesh_boolean." << label << ": "
      << tictoc() << std::endl;
  };
  tictoc();
#endif

  typedef typename DerivedVC::Scalar Scalar;
  //typedef typename DerivedFC::Scalar Index;
  typedef CGAL::Epeck Kernel;
  typedef Kernel::FT ExactScalar;
  typedef Eigen::Matrix<Scalar,Eigen::Dynamic,3> MatrixX3S;
  //typedef Eigen::Matrix<Index,Eigen::Dynamic,Eigen::Dynamic> MatrixXI;
  typedef Eigen::Matrix<typename DerivedJ::Scalar,Eigen::Dynamic,1> VectorXJ;

  // Generate combined mesh.
  typedef Eigen::Matrix<
    ExactScalar,
    Eigen::Dynamic,
    Eigen::Dynamic,
    DerivedVC::IsRowMajor> MatrixXES;
  MatrixXES V;
  DerivedFC F;
  VectorXJ  CJ;
  {
      DerivedVA VV(VA.rows() + VB.rows(), 3);
      DerivedFC FF(FA.rows() + FB.rows(), 3);
      VV << VA, VB;
      FF << FA, FB.array() + VA.rows();
      //// Handle annoying empty cases
      //if(VA.size()>0)
      //{
      //  VV<<VA;
      //}
      //if(VB.size()>0)
      //{
      //  VV<<VB;
      //}
      //if(FA.size()>0)
      //{
      //  FF<<FA;
      //}
      //if(FB.size()>0)
      //{
      //  FF<<FB.array()+VA.rows();
      //}
      resolve_fun(VV, FF, V, F, CJ);
  }
#ifdef MESH_BOOLEAN_TIMING
  log_time("resolve_self_intersection");
#endif

  // Compute winding numbers on each side of each facet.
  const size_t num_faces = F.rows();
  Eigen::MatrixXi W;
  Eigen::VectorXi labels(num_faces);
  std::transform(CJ.data(), CJ.data()+CJ.size(), labels.data(),
      [&](int i) { return i<FA.rows() ? 0:1; });
  bool valid = true;
  if (num_faces > 0) 
  {
    valid = valid & 
      igl::copyleft::cgal::propagate_winding_numbers(V, F, labels, W);
  } else 
  {
    W.resize(0, 4);
  }
  assert((size_t)W.rows() == num_faces);
  if (W.cols() == 2) 
  {
    assert(FB.rows() == 0);
    Eigen::MatrixXi W_tmp(num_faces, 4);
    W_tmp << W, Eigen::MatrixXi::Zero(num_faces, 2);
    W = W_tmp;
  } else {
    assert(W.cols() == 4);
  }
#ifdef MESH_BOOLEAN_TIMING
  log_time("propagate_input_winding_number");
//.........这里部分代码省略.........

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::copyleft::boolean::mesh_boolean(
    const Eigen::PlainObjectBase<DerivedVA> & VA,
    const Eigen::PlainObjectBase<DerivedFA> & FA,
    const Eigen::PlainObjectBase<DerivedVB> & VB,
    const Eigen::PlainObjectBase<DerivedFB> & FB,
    const WindingNumberOp& wind_num_op,
    const KeepFunc& keep,
    const ResolveFunc& resolve_fun,
    Eigen::PlainObjectBase<DerivedVC > & VC,
    Eigen::PlainObjectBase<DerivedFC > & FC,
    Eigen::PlainObjectBase<DerivedJ > & J) {

  typedef typename DerivedVC::Scalar Scalar;
  //typedef typename DerivedFC::Scalar Index;
  typedef CGAL::Epeck Kernel;
  typedef Kernel::FT ExactScalar;
  typedef Eigen::Matrix<Scalar,Eigen::Dynamic,3> MatrixX3S;
  //typedef Eigen::Matrix<Index,Eigen::Dynamic,Eigen::Dynamic> MatrixXI;
  typedef Eigen::Matrix<typename DerivedJ::Scalar,Eigen::Dynamic,1> VectorXJ;

  // Generate combined mesh.
  typedef Eigen::Matrix<
    ExactScalar,
    Eigen::Dynamic,
    Eigen::Dynamic,
    DerivedVC::IsRowMajor> MatrixXES;
  MatrixXES V;
  DerivedFC F;
  VectorXJ  CJ;
  {
      DerivedVA VV(VA.rows() + VB.rows(), 3);
      DerivedFC FF(FA.rows() + FB.rows(), 3);
      VV << VA, VB;
      FF << FA, FB.array() + VA.rows();
      //// Handle annoying empty cases
      //if(VA.size()>0)
      //{
      //  VV<<VA;
      //}
      //if(VB.size()>0)
      //{
      //  VV<<VB;
      //}
      //if(FA.size()>0)
      //{
      //  FF<<FA;
      //}
      //if(FB.size()>0)
      //{
      //  FF<<FB.array()+VA.rows();
      //}
      resolve_fun(VV, FF, V, F, CJ);
  }

  // Compute winding numbers on each side of each facet.
  const size_t num_faces = F.rows();
  Eigen::MatrixXi W;
  Eigen::VectorXi labels(num_faces);
  std::transform(CJ.data(), CJ.data()+CJ.size(), labels.data(),
      [&](int i) { return i<FA.rows() ? 0:1; });
  igl::copyleft::cgal::propagate_winding_numbers(V, F, labels, W);
  assert((size_t)W.rows() == num_faces);
  if (W.cols() == 2) {
    assert(FB.rows() == 0);
    Eigen::MatrixXi W_tmp(num_faces, 4);
    W_tmp << W, Eigen::MatrixXi::Zero(num_faces, 2);
    W = W_tmp;
  } else {
    assert(W.cols() == 4);
  }

  // Compute resulting winding number.
  Eigen::MatrixXi Wr(num_faces, 2);
  for (size_t i=0; i<num_faces; i++) {
    Eigen::MatrixXi w_out(1,2), w_in(1,2);
    w_out << W(i,0), W(i,2);
    w_in  << W(i,1), W(i,3);
    Wr(i,0) = wind_num_op(w_out);
    Wr(i,1) = wind_num_op(w_in);
  }

  // Extract boundary separating inside from outside.
  auto index_to_signed_index = [&](size_t i, bool ori) -> int{
    return (i+1)*(ori?1:-1);
  };
  //auto signed_index_to_index = [&](int i) -> size_t {
  //    return abs(i) - 1;
  //};
  std::vector<int> selected;
  for(size_t i=0; i<num_faces; i++) {
    auto should_keep = keep(Wr(i,0), Wr(i,1));
    if (should_keep > 0) {
      selected.push_back(index_to_signed_index(i, true));
    } else if (should_keep < 0) {
      selected.push_back(index_to_signed_index(i, false));
    }
  }

  const size_t num_selected = selected.size();
  DerivedFC kept_faces(num_selected, 3);
//.........这里部分代码省略.........