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

IGL_INLINE void igl::sortrows(
  const Eigen::PlainObjectBase<DerivedX>& X,
  const bool ascending,
  Eigen::PlainObjectBase<DerivedX>& Y,
  Eigen::PlainObjectBase<DerivedIX>& IX)
{
  using namespace std;
  using namespace Eigen;
  using namespace igl;
  // Resize output
  Y.resize(X.rows(),X.cols());
  IX.resize(X.rows(),1);
  for(int i = 0;i<X.rows();i++)
  {
    IX(i) = i;
  }
  std::sort(
    IX.data(),
    IX.data()+IX.size(),
    igl::IndexRowLessThan<const Eigen::PlainObjectBase<DerivedX> & >(X));
  // if not ascending then reverse
  if(!ascending)
  {
    std::reverse(IX.data(),IX.data()+IX.size());
  }
  for(int i = 0;i<X.rows();i++)
  {
    Y.row(i) = X.row(IX(i));
  }
}

IGL_INLINE void igl::setdiff(
  const Eigen::PlainObjectBase<DerivedA> & A,
  const Eigen::PlainObjectBase<DerivedB> & B,
  Eigen::PlainObjectBase<DerivedC> & C,
  Eigen::PlainObjectBase<DerivedIA> & IA)
{
  using namespace Eigen;
  using namespace std;
  // boring base cases
  if(A.size() == 0)
  {
    C.resize(0,1);
    IA.resize(0,1);
  }
  if(B.size() == 0)
  {
    C.resize(A.size(),1);
    copy(A.data(),A.data()+A.size(),C.data());
    IA = igl::colon<typename DerivedIA::Scalar>(0,C.size()-1);
  }

  // Get rid of any duplicates
  typedef Matrix<typename DerivedA::Scalar,Dynamic,1> VectorA;
  typedef Matrix<typename DerivedB::Scalar,Dynamic,1> VectorB;
  VectorA uA;
  VectorB uB;
  typedef PlainObjectBase<DerivedIA> IAType;
  IAType uIA,uIuA,uIB,uIuB;
  unique(A,uA,uIA,uIuA);
  unique(B,uB,uIB,uIuB);

  // Sort both
  VectorA sA;
  VectorB sB;
  IAType sIA,sIB;
  sort(uA,1,true,sA,sIA);
  sort(uB,1,true,sB,sIB);

  vector<typename DerivedB::Scalar> vC;
  vector<typename DerivedIA::Scalar> vIA;
  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))
    {
      // then sA(a) did not appear in sB
      vC.push_back(sA(a));
      vIA.push_back(uIA(sIA(a)));
    }
  }
  list_to_matrix(vC,C);
  list_to_matrix(vIA,IA);
}

IGL_INLINE void igl::randperm(
  const int n,
  Eigen::PlainObjectBase<DerivedI> & I)
{
  Eigen::VectorXi II;
  igl::colon(0,1,n-1,II);
  I = II;
  std::random_shuffle(I.data(),I.data()+n);
}

IGL_INLINE int igl::unproject_to_zero_plane(
  const Eigen::PlainObjectBase<Derivedwin> & win,
  Eigen::PlainObjectBase<Derivedobj> & obj)
{
  return unproject_to_zero_plane(win(0),win(1),
      &obj.data()[0],
      &obj.data()[1],
      &obj.data()[2]);
}

void clamp(Eigen::PlainObjectBase<Derived>& m, typename Derived::Scalar lower, typename Derived::Scalar upper) {
    typedef typename Eigen::PlainObjectBase<Derived>::Index idx_t;

    for (idx_t i = 0; i < m.outerSize(); ++i) {
        for (idx_t j = 0; j < m.innerSize(); ++j) {
            int offset = i * m.innerSize() + j;
            if (m.data()[offset] < lower) m.data()[offset] = lower;
            if (m.data()[offset] > upper) m.data()[offset] = upper;
        }
    }
}

IGL_INLINE void igl::mod(
  const Eigen::PlainObjectBase<DerivedA> & A,
  const int base,
  Eigen::PlainObjectBase<DerivedB> & B)
{
  B.resize(A.rows(),A.cols());
  for(int i = 0;i<B.size();i++)
  {
    *(B.data()+i) = (*(A.data()+i))%base;
  }
}

IGL_INLINE void igl::opengl::vertex_array(
  const Eigen::PlainObjectBase<DerivedV> & V,
  const Eigen::PlainObjectBase<DerivedF> & F,
  GLuint & va_id,
  GLuint & ab_id,
  GLuint & eab_id)
{
  // Inputs should be in RowMajor storage. If not, we have no choice but to
  // create a copy.
  if(!(V.Options & Eigen::RowMajor))
  {
    Eigen::Matrix<
      typename DerivedV::Scalar,
      DerivedV::RowsAtCompileTime,
      DerivedV::ColsAtCompileTime,
      Eigen::RowMajor> VR = V;
    return vertex_array(VR,F,va_id,ab_id,eab_id);
  }
  if(!(F.Options & Eigen::RowMajor))
  {
    Eigen::Matrix<
      typename DerivedF::Scalar,
      DerivedF::RowsAtCompileTime,
      DerivedF::ColsAtCompileTime,
      Eigen::RowMajor> FR = F;
    return vertex_array(V,FR,va_id,ab_id,eab_id);
  }
  // Generate and attach buffers to vertex array
  glGenVertexArrays(1, &va_id);
  glGenBuffers(1, &ab_id);
  glGenBuffers(1, &eab_id);
  glBindVertexArray(va_id);
  glBindBuffer(GL_ARRAY_BUFFER, ab_id);
  const auto size_VScalar = sizeof(typename DerivedV::Scalar);
  const auto size_FScalar = sizeof(typename DerivedF::Scalar);
  glBufferData(GL_ARRAY_BUFFER,size_VScalar*V.size(),V.data(),GL_STATIC_DRAW);
  glBindBuffer(GL_ELEMENT_ARRAY_BUFFER, eab_id);
  assert(sizeof(GLuint) == size_FScalar && "F type does not match GLuint");
  glBufferData(
    GL_ELEMENT_ARRAY_BUFFER, sizeof(GLuint)*F.size(), F.data(), GL_STATIC_DRAW);
  glVertexAttribPointer(
    0,
    V.cols(),
    size_VScalar==sizeof(float)?GL_FLOAT:GL_DOUBLE,
    GL_FALSE,
    V.cols()*size_VScalar,
    (GLvoid*)0);
  glEnableVertexAttribArray(0);
  glBindBuffer(GL_ARRAY_BUFFER, 0); 
  glBindVertexArray(0);
}

IGL_INLINE void igl::matlab::parse_rhs_double(
    const mxArray *prhs[], 
    Eigen::PlainObjectBase<DerivedV> & V)
{
  using namespace std;
  using namespace Eigen;
  // set number of mesh vertices
  const int n = mxGetM(prhs[0]);
  // set vertex position pointers
  double * Vp = mxGetPr(prhs[0]);
  const int dim = mxGetN(prhs[0]);

  typedef typename DerivedV::Scalar Scalar;
  Matrix<Scalar, DerivedV::ColsAtCompileTime, DerivedV::RowsAtCompileTime, RowMajor> VT;
  Scalar * V_data;
  if(DerivedV::IsRowMajor)
  {
    VT.resize(dim,n);
    V_data = VT.data();
  }else
  {
    V.resize(n,dim);
    V_data = V.data();
  }
  copy(Vp,Vp+n*dim,V_data);
  if(DerivedV::IsRowMajor)
  {
    V = VT.transpose();
  }
}

IGL_INLINE void igl::outer_vertex(
        const Eigen::PlainObjectBase<DerivedV> & V,
        const Eigen::PlainObjectBase<DerivedF> & F,
        const Eigen::PlainObjectBase<DerivedI> & I,
        IndexType & v_index,
        Eigen::PlainObjectBase<DerivedA> & A) {
    // Algorithm: 
    //    Find an outer vertex (i.e. vertex reachable from infinity)
    //    Return the vertex with the largest X value.
    //    If there is a tie, pick the one with largest Y value.
    //    If there is still a tie, pick the one with the largest Z value.
    //    If there is still a tie, then there are duplicated vertices within the
    //    mesh, which violates the precondition.
    const size_t INVALID = std::numeric_limits<size_t>::max();
    const size_t num_selected_faces = I.rows();
    std::vector<size_t> candidate_faces;
    size_t outer_vid = INVALID;
    typename DerivedV::Scalar outer_val = 0;
    for (size_t i=0; i<num_selected_faces; i++) {
        size_t f = I(i);
        for (size_t j=0; j<3; j++) {
            auto v = F(f, j);
            auto vx = V(v, 0);
            if (outer_vid == INVALID || vx > outer_val) {
                outer_val = vx;
                outer_vid = v;
                candidate_faces = {f};
            } else if (v == outer_vid) {
                candidate_faces.push_back(f);
            } else if (vx == outer_val) {
                // Break tie.
                auto vy = V(v,1);
                auto vz = V(v, 2);
                auto outer_y = V(outer_vid, 1);
                auto outer_z = V(outer_vid, 2);
                assert(!(vy == outer_y && vz == outer_z));
                bool replace = (vy > outer_y) ||
                    ((vy == outer_y) && (vz > outer_z));
                if (replace) {
                    outer_val = vx;
                    outer_vid = v;
                    candidate_faces = {f};
                }
            }
        }
    }

    assert(outer_vid != INVALID);
    assert(candidate_faces.size() > 0);
    v_index = outer_vid;
    A.resize(candidate_faces.size());
    std::copy(candidate_faces.begin(), candidate_faces.end(), A.data());
}

IGL_INLINE void igl::parse_rhs_double(
    const mxArray *prhs[], 
    Eigen::PlainObjectBase<DerivedV> & V)
{
  using namespace std;
  // set number of mesh vertices
  const int n = mxGetM(prhs[0]);
  // set vertex position pointers
  double * Vp = mxGetPr(prhs[0]);
  const int dim = mxGetN(prhs[0]);
  V.resize(n,dim);
  copy(Vp,Vp+n*dim,&V.data()[0]);
}

IGL_INLINE void igl::components(
  const Eigen::SparseMatrix<AScalar> & A,
  Eigen::PlainObjectBase<DerivedC> & C)
{
  assert(A.rows() == A.cols());
  using namespace Eigen;
  // THIS IS DENSE:
  //boost::adjacency_matrix<boost::undirectedS> bA(A.rows());
  boost::adjacency_list<boost::vecS,boost::vecS,boost::undirectedS> bA(A.rows());
  for(int j=0; j<A.outerSize();j++)
  {
    // Iterate over inside
    for(typename SparseMatrix<AScalar>::InnerIterator it (A,j); it; ++it)
    {
      if(0 != it.value())
      {
        boost::add_edge(it.row(),it.col(),bA);
      }
    }
  }
  C.resize(A.rows(),1);
  boost::connected_components(bA,C.data());
}

IGL_INLINE void igl::resolve_duplicated_faces(
    const Eigen::PlainObjectBase<DerivedF1>& F1,
    Eigen::PlainObjectBase<DerivedF2>& F2,
    Eigen::PlainObjectBase<DerivedJ>& J) {

  //typedef typename DerivedF1::Scalar Index;
  Eigen::VectorXi IA,IC;
  DerivedF1 uF;
  igl::unique_simplices(F1,uF,IA,IC);

  const size_t num_faces = F1.rows();
  const size_t num_unique_faces = uF.rows();
  assert((size_t) IA.rows() == num_unique_faces);
  // faces on top of each unique face
  std::vector<std::vector<int> > uF2F(num_unique_faces);
  // signed counts
  Eigen::VectorXi counts = Eigen::VectorXi::Zero(num_unique_faces);
  Eigen::VectorXi ucounts = Eigen::VectorXi::Zero(num_unique_faces);
  // loop over all faces
  for (size_t i=0; i<num_faces; i++) {
    const size_t ui = IC(i);
    const bool consistent = 
      (F1(i,0) == uF(ui, 0) && F1(i,1) == uF(ui, 1) && F1(i,2) == uF(ui, 2)) ||
      (F1(i,0) == uF(ui, 1) && F1(i,1) == uF(ui, 2) && F1(i,2) == uF(ui, 0)) ||
      (F1(i,0) == uF(ui, 2) && F1(i,1) == uF(ui, 0) && F1(i,2) == uF(ui, 1));
    uF2F[ui].push_back(int(i+1) * (consistent?1:-1));
    counts(ui) += consistent ? 1:-1;
    ucounts(ui)++;
  }

  std::vector<size_t> kept_faces;
  for (size_t i=0; i<num_unique_faces; i++) {
    if (ucounts[i] == 1) {
      kept_faces.push_back(abs(uF2F[i][0])-1);
      continue;
    }
    if (counts[i] == 1) {
      bool found = false;
      for (auto fid : uF2F[i]) {
        if (fid > 0) {
          kept_faces.push_back(abs(fid)-1);
          found = true;
          break;
        }
      }
      assert(found);
    } else if (counts[i] == -1) {
      bool found = false;
      for (auto fid : uF2F[i]) {
        if (fid < 0) {
          kept_faces.push_back(abs(fid)-1);
          found = true;
          break;
        }
      }
      assert(found);
    } else {
      assert(counts[i] == 0);
    }
  }

  const size_t num_kept = kept_faces.size();
  J.resize(num_kept, 1);
  std::copy(kept_faces.begin(), kept_faces.end(), J.data());
  igl::slice(F1, J, 1, F2);
}

//.........这里部分代码省略.........
    {
      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();
      }
    }
  }
}

IGL_INLINE void igl::copyleft::cgal::snap_rounding(
  const Eigen::PlainObjectBase<DerivedV> & V,
  const Eigen::PlainObjectBase<DerivedE> & E,
  Eigen::PlainObjectBase<DerivedVI> & VI,
  Eigen::PlainObjectBase<DerivedEI> & EI,
  Eigen::PlainObjectBase<DerivedJ> & J)
{
  using namespace Eigen;
  using namespace igl;
  using namespace igl::copyleft::cgal;
  using namespace std;
  // Exact scalar type
  typedef CGAL::Epeck Kernel;
  typedef Kernel::FT EScalar;
  typedef CGAL::Segment_2<Kernel> Segment_2;
  typedef CGAL::Point_2<Kernel> Point_2;
  typedef CGAL::Vector_2<Kernel> Vector_2;
  typedef Matrix<EScalar,Dynamic,Dynamic>  MatrixXE;
  // Convert vertex positions to exact kernel

  MatrixXE VE;
  {
    VectorXi IM;
    resolve_intersections(V,E,VE,EI,J,IM);
    for_each(
      EI.data(),
      EI.data()+EI.size(),
      [&IM](typename DerivedEI::Scalar& i){i=IM(i);});
    VectorXi _;
    remove_unreferenced( MatrixXE(VE), DerivedEI(EI), VE,EI,_);
  }

  // find all hot pixels
  //// southwest and north east corners
  //const RowVector2i SW(
  //  round(VE.col(0).minCoeff()),
  //  round(VE.col(1).minCoeff()));
  //const RowVector2i NE(
  //  round(VE.col(0).maxCoeff()),
  //  round(VE.col(1).maxCoeff()));

  // https://github.com/CGAL/cgal/issues/548
  // Round an exact scalar to the nearest integer. A priori can't just round
  // double. Suppose e=0.5+ε but double(e)<0.5
  //
  // Inputs:
  //   e  exact number
  // Outputs:
  //   i  closest integer
  const auto & round = [](const EScalar & e)->int
  {
    const double d = CGAL::to_double(e);
    // get an integer that's near the closest int
    int i = std::round(d);
    EScalar di_sqr = CGAL::square((e-EScalar(i)));
    const auto & search = [&i,&di_sqr,&e](const int dir)
    {
      while(true)
      {
        const int j = i+dir;
        const EScalar dj_sqr = CGAL::square((e-EScalar(j)));
        if(dj_sqr < di_sqr)
        {
          i = j;
          di_sqr = dj_sqr;
        }else
        {
          break;
        }
      }
    };
    // Try to increase/decrease int
    search(1);
    search(-1);
    return i;
  };
  vector<Point_2> hot;
  for(int i = 0;i<VE.rows();i++)
  {
    hot.emplace_back(round(VE(i,0)),round(VE(i,1)));
  }
  {
    std::vector<size_t> _1,_2;
    igl::unique(vector<Point_2>(hot),hot,_1,_2);
  }

  // find all segments intersecting hot pixels
  //   split edge at closest point to hot pixel center
  vector<vector<Point_2>>  steiner(EI.rows());
  // initialize each segment with endpoints
  for(int i = 0;i<EI.rows();i++)
  {
    steiner[i].emplace_back(VE(EI(i,0),0),VE(EI(i,0),1));
    steiner[i].emplace_back(VE(EI(i,1),0),VE(EI(i,1),1));
  }
  // silly O(n²) implementation
  for(const Point_2 & h : hot)
  {
    // North, East, South, West
    Segment_2 wall[4] = 
//.........这里部分代码省略.........

IGL_INLINE void igl::outer_edge(
        const Eigen::PlainObjectBase<DerivedV> & V,
        const Eigen::PlainObjectBase<DerivedF> & F,
        const Eigen::PlainObjectBase<DerivedI> & I,
        IndexType & v1,
        IndexType & v2,
        Eigen::PlainObjectBase<DerivedA> & A) {
    // Algorithm:
    //    Find an outer vertex first.
    //    Find the incident edge with largest slope when projected onto XY plane.
    //    If there is still a tie, break it using the projected slope onto ZX plane.
    //    If there is still a tie, then there are multiple overlapping edges,
    //    which violates the precondition.
    typedef typename DerivedV::Scalar Scalar;
    typedef typename DerivedV::Index Index;
    typedef typename Eigen::Matrix<Scalar, 3, 1> ScalarArray3;
    typedef typename Eigen::Matrix<typename DerivedF::Scalar, 3, 1> IndexArray3;
    const size_t INVALID = std::numeric_limits<size_t>::max();

    Index outer_vid;
    Eigen::Matrix<Index,Eigen::Dynamic,1> candidate_faces;
    outer_vertex(V, F, I, outer_vid, candidate_faces);
    const ScalarArray3& outer_v = V.row(outer_vid);
    assert(candidate_faces.size() > 0);

    auto get_vertex_index = [&](const IndexArray3& f, Index vid) -> Index 
    {
        if (f[0] == vid) return 0;
        if (f[1] == vid) return 1;
        if (f[2] == vid) return 2;
        assert(false);
        return -1;
    };

    Scalar outer_slope_YX = 0;
    Scalar outer_slope_ZX = 0;
    size_t outer_opp_vid = INVALID;
    bool infinite_slope_detected = false;
    std::vector<Index> incident_faces;
    auto check_and_update_outer_edge = [&](Index opp_vid, Index fid) {
        if (opp_vid == outer_opp_vid) {
            incident_faces.push_back(fid);
            return;
        }

        const ScalarArray3 opp_v = V.row(opp_vid);
        if (!infinite_slope_detected && outer_v[0] != opp_v[0]) {
            // Finite slope
            const ScalarArray3 diff = opp_v - outer_v;
            const Scalar slope_YX = diff[1] / diff[0];
            const Scalar slope_ZX = diff[2] / diff[0];
            if (outer_opp_vid == INVALID ||
                    slope_YX > outer_slope_YX ||
                    (slope_YX == outer_slope_YX &&
                     slope_ZX > outer_slope_ZX)) {
                outer_opp_vid = opp_vid;
                outer_slope_YX = slope_YX;
                outer_slope_ZX = slope_ZX;
                incident_faces = {fid};
            }
        } else if (!infinite_slope_detected) {
            // Infinite slope
            outer_opp_vid = opp_vid;
            infinite_slope_detected = true;
            incident_faces = {fid};
        }
    };

    const auto num_candidate_faces = candidate_faces.size();
    for (size_t i=0; i<num_candidate_faces; i++) {
        const Index fid = candidate_faces(i);
        const IndexArray3& f = F.row(fid);
        size_t id = get_vertex_index(f, outer_vid);
        Index next_vid = f((id+1)%3);
        Index prev_vid = f((id+2)%3);
        check_and_update_outer_edge(next_vid, fid);
        check_and_update_outer_edge(prev_vid, fid);
    }

    v1 = outer_vid;
    v2 = outer_opp_vid;
    A.resize(incident_faces.size());
    std::copy(incident_faces.begin(), incident_faces.end(), A.data());
}