C++ MatrixXi::row方法代码示例

void relabelFaces(Eigen::MatrixXi& aggregated,
                  const Eigen::MatrixXd& vertices,
                  const Eigen::MatrixXi& faces,
                  Eigen::Vector3i& vertexOffset,
                  int& offset) {
  for (int i = 0; i < faces.rows(); i++) {
    aggregated.row(offset++) = vertexOffset + Eigen::Vector3i(faces.row(i));
  }
  int numVerts = vertices.rows();
  vertexOffset += Eigen::Vector3i(numVerts, numVerts, numVerts);
}

bool Polygon::convertToInequalityConstraints(Eigen::MatrixXd& A, Eigen::VectorXd& b) const
{
  Eigen::MatrixXd V(nVertices(), 2);
  for (unsigned int i = 0; i < nVertices(); ++i)
    V.row(i) = vertices_[i];

  // Create k, a list of indices from V forming the convex hull.
  // TODO: Assuming counter-clockwise ordered convex polygon.
  // MATLAB: k = convhulln(V);
  Eigen::MatrixXi k;
  k.resizeLike(V);
  for (unsigned int i = 0; i < V.rows(); ++i)
    k.row(i) << i, (i+1) % V.rows();
  Eigen::RowVectorXd c = V.colwise().mean();
  V.rowwise() -= c;
  A = Eigen::MatrixXd::Constant(k.rows(), V.cols(), NAN);

  unsigned int rc = 0;
  for (unsigned int ix = 0; ix < k.rows(); ++ix) {
    Eigen::MatrixXd F(2, V.cols());
    F.row(0) << V.row(k(ix, 0));
    F.row(1) << V.row(k(ix, 1));
    Eigen::FullPivLU<Eigen::MatrixXd> luDecomp(F);
    if (luDecomp.rank() == F.rows()) {
      A.row(rc) = F.colPivHouseholderQr().solve(Eigen::VectorXd::Ones(F.rows()));
      ++rc;
    }
  }

  A = A.topRows(rc);
  b = Eigen::VectorXd::Ones(A.rows());
  b = b + A * c.transpose();

  return true;
}

void shortest_edge_and_midpoint1(const int e, const Eigen::MatrixXd &V,
                                 const Eigen::MatrixXi &F,
                                 const Eigen::MatrixXi &E,
                                 const Eigen::VectorXi &EMAP,
                                 const Eigen::MatrixXi &EF,
                                 const Eigen::MatrixXi &EI, double &cost,
                                 RowVectorXd &p) {
    // vectorsum
    const int eflip = E(e, 0) > E(e, 1);
    const std::vector<int> nV2Fd =
        igl::circulation(e, !eflip, F, E, EMAP, EF, EI);
    p = 0.5 * (V.row(E(e, 0)) + V.row(E(e, 1)));
    Eigen::RowVectorXd pointy(3);
    pointy.setZero();
    std::set<int> newEdges;
    for (int i = 0; i < nV2Fd.size(); i++) {
        for (int j = 0; j < 3; j++) {
            int curVert = F.row(nV2Fd[i])[j];
            if (curVert != E(e, 0) || curVert != E(e, 1)) {
                if (newEdges.insert(curVert).second) {
                    pointy = (V.row(curVert) - p) + pointy;
                }
            }
        }
    }
    cost = (pointy).norm();
}

int main(int argc, char *argv[])
{
  using namespace Eigen;
  using namespace std;

  cout<<"Usage:"<<endl;
  cout<<"[space]  toggle showing input mesh, output mesh or slice "<<endl;
  cout<<"         through tet-mesh of convex hull."<<endl;
  cout<<"'.'/','  push back/pull forward slicing plane."<<endl;
  cout<<endl;

  // Load mesh: (V,T) tet-mesh of convex hull, F contains facets of input
  // surface mesh _after_ self-intersection resolution
  igl::readMESH(TUTORIAL_SHARED_PATH "/big-sigcat.mesh",V,T,F);

  // Compute barycenters of all tets
  igl::barycenter(V,T,BC);

  // Compute generalized winding number at all barycenters
  cout<<"Computing winding number over all "<<T.rows()<<" tets..."<<endl;
  igl::winding_number(V,F,BC,W);

  // Extract interior tets
  MatrixXi CT((W.array()>0.5).count(),4);
  {
    size_t k = 0;
    for(size_t t = 0;t<T.rows();t++)
    {
      if(W(t)>0.5)
      {
        CT.row(k) = T.row(t);
        k++;
      }
    }
  }
  // find bounary facets of interior tets
  igl::boundary_facets(CT,G);
  // boundary_facets seems to be reversed...
  G = G.rowwise().reverse().eval();

  // normalize
  W = (W.array() - W.minCoeff())/(W.maxCoeff()-W.minCoeff());

  // Plot the generated mesh
  igl::opengl::glfw::Viewer viewer;
  update_visualization(viewer);
  viewer.callback_key_down = &key_down;
  viewer.launch();
}

void igl::polyvector_field_one_ring_matchings(const Eigen::PlainObjectBase<DerivedV> &V,
                                              const Eigen::PlainObjectBase<DerivedF> &F,
                                              const std::vector<std::vector<VFType> >& VF,
                                              const Eigen::MatrixXi& E2F,
                                              const Eigen::MatrixXi& F2E,
                                              const Eigen::PlainObjectBase<DerivedTT>& TT,
                                              const Eigen::PlainObjectBase<DerivedM> &match_ab,
                                              const Eigen::PlainObjectBase<DerivedM> &match_ba,
                                              const int vi,
                                              Eigen::MatrixXi &mvi,
                                              Eigen::VectorXi &fi)
{
  int half_degree = match_ab.cols();
  mvi.resize(VF[vi].size()+1,half_degree);
  fi.resize(VF[vi].size()+1,1);
  //start from one face
  //first, check if the vertex is on a boundary
  //then there must be two faces that are on the boundary
  //(other cases not supported)

  int fstart = -1;
  int ind = 0;
  for (int i =0; i<VF[vi].size(); ++i)
  {
    int fi = VF[vi][i];
    for (int  j=0; j<3; ++j)
      if (F(fi,j)==vi && TT(fi,j) == -1)
      {
        ind ++;
        fstart = fi;
        //        break;
      }
  }
  if (ind >1 )
  {
    std::cerr<<"igl::polyvector_field_one_ring_matchings -- vertex "<<vi<< " is on an unusual boundary"<<std::endl;
    exit(1);
  }
  if (fstart == -1)
    fstart = VF[vi][0];
  int current_face = fstart;
  int i =0;
  fi[i] = current_face;
  for (int j=0; j<half_degree; ++j)
    mvi(i,j) = j;

  int next_face = -1;
  while (next_face != fstart && current_face!=-1)
  {
    // look for the vertex
    int j=-1;
    for (unsigned z=0; z<3; ++z)
      if (F(current_face,(z+1)%3) == vi)
      {
        j=z;
        break;
      }
    assert(j!=-1);

    next_face = TT(current_face, j);
    ++i;

    if (next_face == -1)
      mvi.row(i).setConstant(-1);
    else
    {
      // look at the edge between the two faces
      const int &current_edge = F2E(current_face,j);

      for (int k=0; k<half_degree; ++k)
      {
        // check its orientation to determine whether match_ab or match_ba should be used
        if ((E2F(current_edge,0) == current_face) &&
            (E2F(current_edge,1) == next_face) )
        {
          //look at match_ab
          mvi(i,k) = match_ab(current_edge,(mvi(i-1,k))%half_degree);
        }
        else
        {
          assert((E2F(current_edge,1) == current_face) &&
                 (E2F(current_edge,0) == next_face));
          //look at match_ba
          mvi(i,k) = match_ba(current_edge,(mvi(i-1,k))%half_degree);
        }
        if (mvi(i-1,k)>=half_degree)
          mvi(i,k) = (mvi(i,k)+half_degree)%(2*half_degree);
      }
    }
    current_face = next_face;
    fi[i] = current_face;
  }
}

void igl::collapse_small_triangles(
  const Eigen::MatrixXd & V,
  const Eigen::MatrixXi & F,
  const double eps,
  Eigen::MatrixXi & FF)
{
  using namespace Eigen;
  using namespace std;

  // Compute bounding box diagonal length
  double bbd = bounding_box_diagonal(V);
  MatrixXd l;
  edge_lengths(V,F,l);
  VectorXd dblA;
  doublearea(l,dblA);

  // Minimum area tolerance
  const double min_dblarea = 2.0*eps*bbd*bbd;

  Eigen::VectorXi FIM = colon<int>(0,V.rows()-1);
  int num_edge_collapses = 0;
  // Loop over triangles
  for(int f = 0;f<F.rows();f++)
  {
    if(dblA(f) < min_dblarea)
    {
      double minl = 0;
      int minli = -1;
      // Find shortest edge
      for(int e = 0;e<3;e++)
      {
        if(minli==-1 || l(f,e)<minl)
        {
          minli = e;
          minl = l(f,e);
        }
      }
      double maxl = 0;
      int maxli = -1;
      // Find longest edge
      for(int e = 0;e<3;e++)
      {
        if(maxli==-1 || l(f,e)>maxl)
        {
          maxli = e;
          maxl = l(f,e);
        }
      }
      // Be sure that min and max aren't the same
      maxli = (minli==maxli?(minli+1)%3:maxli);

      // Collapse min edge maintaining max edge: i-->j
      // Q: Why this direction?
      int i = maxli;
      int j = ((minli+1)%3 == maxli ? (minli+2)%3: (minli+1)%3);
      assert(i != minli);
      assert(j != minli);
      assert(i != j);
      FIM(F(f,i)) = FIM(F(f,j));
      num_edge_collapses++;
    }
  }

  // Reindex faces
  MatrixXi rF = F;
  // Loop over triangles
  for(int f = 0;f<rF.rows();f++)
  {
    for(int i = 0;i<rF.cols();i++)
    {
      rF(f,i) = FIM(rF(f,i));
    }
  }

  FF.resize(rF.rows(),rF.cols());
  int num_face_collapses=0;
  // Only keep uncollapsed faces
  {
    int ff = 0;
    // Loop over triangles
    for(int f = 0;f<rF.rows();f++)
    {
      bool collapsed = false;
      // Check if any indices are the same
      for(int i = 0;i<rF.cols();i++)
      {
        for(int j = i+1;j<rF.cols();j++)
        {
          if(rF(f,i)==rF(f,j))
          {
            collapsed = true;
            num_face_collapses++;
            break;
          }
        }
      }
      if(!collapsed)
      {
        FF.row(ff++) = rF.row(f);
      }
//.........这里部分代码省略.........

// A tet is organized in the following fashion:
//
//       0
//       
//       3       (3 is in the background; 1 and 2 are in the foreground)
//  2        1
//
// So the faces (with counterclockwise normals) are:
// (0 1 3)
// (0 2 1)
// (3 2 0)
// (1 2 3)
//
//
// This method will perform three steps:
// 1. Add all faces, duplicating ones on the interior
// 2. Remove all duplicate verts (might have been caused during #1)
// 3. Remove all duplicate faces
void marching_tets(
  const Eigen::MatrixXd& V,
  const Eigen::MatrixXi& T,
  const Eigen::VectorXd& H,
  double offset,
  Eigen::MatrixXd& NV,
  Eigen::MatrixXi& NF,
  Eigen::VectorXi& I)
{
  min_how_much = 1;
  max_how_much = 0;
  using namespace Eigen;
  using namespace std;

  // Count the faces.
  std::map<std::vector<int>, int> face_counts;
  for (int i = 0; i < T.rows(); ++i) {
    std::vector<std::vector<int> > fs;
    fs.push_back({T(i, 0), T(i, 1), T(i, 3)});
    fs.push_back({T(i, 0), T(i, 2), T(i, 1)});
    fs.push_back({T(i, 3), T(i, 2), T(i, 0)});
    fs.push_back({T(i, 1), T(i, 2), T(i, 3)});
    for (auto &f : fs) {
      std::sort(f.begin(), f.end());
      // Add it to the map.
      face_counts[f]++;
    }
  }

  vector<Eigen::RowVector3i> faces;
  vector<int> faces_markers;
  vector<Eigen::RowVector3d> new_verts;
  int times[6];
  for (int i = 0; i < 6; i++) {
    times[i] = 0;
  }

  // Create data structure.
  MarchingTetsDat dd(V, H, faces, faces_markers, face_counts, new_verts, offset);

  int numEq = 0;

  // Check each tet face, add as needed.
  for (int i = 0; i < T.rows(); ++i) {
    // See if the tet is entirely inside.
    vector<int> inside, outside, inside_t, outside_t, identical;
    for (int j = 0; j < T.cols(); ++j) {
      //if (H(T(i, j)) > offset + 1e-4) {
      if (H(T(i, j)) > offset) {
        outside.push_back(j);
      } else if (H(T(i, j)) < offset) {
        inside.push_back(j);
      } else {
        numEq++;
        identical.push_back(j);
      }

      if (H(T(i, j)) == GLOBAL::outside_temp) {
        outside_t.push_back(j);
      } else if (H(T(i, j)) == GLOBAL::inside_temp) {
        inside_t.push_back(j);
      }
    }

    // Ignore this tet if it's entirely outside.
    if (outside.size() == 4) {
      continue;
    }


    if (outside.size() == 0 && inside.size() == 0) {
      // degenerate, ignore.
      printf("WARNING: degenerate tet face found!!\n");
    } else if (inside.size() == 0 && identical.size() < 3) {
      // Nothing to add.
    } else if (identical.size() == 3) {
      //addOrig(T.row(i), 7, dd.faces, dd.faces_markers);
      // Ignore it if there's only one on the outside.
      //if (inside.size() == 0) continue;
      if (outside.size() == 0) continue;
      times[1]++;
      // Add just a single face (always)
//.........这里部分代码省略.........

IGL_INLINE bool igl::decimate(
  const Eigen::MatrixXd & OV,
  const Eigen::MatrixXi & OF,
  const std::function<void(
    const int,
    const Eigen::MatrixXd &,
    const Eigen::MatrixXi &,
    const Eigen::MatrixXi &,
    const Eigen::VectorXi &,
    const Eigen::MatrixXi &,
    const Eigen::MatrixXi &,
    double &,
    Eigen::RowVectorXd &)> & cost_and_placement,
  const std::function<bool(
      const Eigen::MatrixXd &,
      const Eigen::MatrixXi &,
      const Eigen::MatrixXi &,
      const Eigen::VectorXi &,
      const Eigen::MatrixXi &,
      const Eigen::MatrixXi &,
      const std::set<std::pair<double,int> > &,
      const std::vector<std::set<std::pair<double,int> >::iterator > &,
      const Eigen::MatrixXd &,
      const int,
      const int,
      const int,
      const int,
      const int)> & stopping_condition,
    const std::function<bool(
      const Eigen::MatrixXd &                                         ,/*V*/
      const Eigen::MatrixXi &                                         ,/*F*/
      const Eigen::MatrixXi &                                         ,/*E*/
      const Eigen::VectorXi &                                         ,/*EMAP*/
      const Eigen::MatrixXi &                                         ,/*EF*/
      const Eigen::MatrixXi &                                         ,/*EI*/
      const std::set<std::pair<double,int> > &                        ,/*Q*/
      const std::vector<std::set<std::pair<double,int> >::iterator > &,/*Qit*/
      const Eigen::MatrixXd &                                         ,/*C*/
      const int                                                        /*e*/
      )> & pre_collapse,
    const std::function<void(
      const Eigen::MatrixXd &                                         ,   /*V*/
      const Eigen::MatrixXi &                                         ,   /*F*/
      const Eigen::MatrixXi &                                         ,   /*E*/
      const Eigen::VectorXi &                                         ,/*EMAP*/
      const Eigen::MatrixXi &                                         ,  /*EF*/
      const Eigen::MatrixXi &                                         ,  /*EI*/
      const std::set<std::pair<double,int> > &                        ,   /*Q*/
      const std::vector<std::set<std::pair<double,int> >::iterator > &, /*Qit*/
      const Eigen::MatrixXd &                                         ,   /*C*/
      const int                                                       ,   /*e*/
      const int                                                       ,  /*e1*/
      const int                                                       ,  /*e2*/
      const int                                                       ,  /*f1*/
      const int                                                       ,  /*f2*/
      const bool                                                  /*collapsed*/
      )> & post_collapse,
  const Eigen::MatrixXi & OE,
  const Eigen::VectorXi & OEMAP,
  const Eigen::MatrixXi & OEF,
  const Eigen::MatrixXi & OEI,
  Eigen::MatrixXd & U,
  Eigen::MatrixXi & G,
  Eigen::VectorXi & J,
  Eigen::VectorXi & I
  )
{

  // Decimate 1
  using namespace Eigen;
  using namespace std;
  // Working copies
  Eigen::MatrixXd V = OV;
  Eigen::MatrixXi F = OF;
  Eigen::MatrixXi E = OE;
  Eigen::VectorXi EMAP = OEMAP;
  Eigen::MatrixXi EF = OEF;
  Eigen::MatrixXi EI = OEI;
  typedef std::set<std::pair<double,int> > PriorityQueue;
  PriorityQueue Q;
  std::vector<PriorityQueue::iterator > Qit;
  Qit.resize(E.rows());
  // If an edge were collapsed, we'd collapse it to these points:
  MatrixXd C(E.rows(),V.cols());
  for(int e = 0;e<E.rows();e++)
  {
    double cost = e;
    RowVectorXd p(1,3);
    cost_and_placement(e,V,F,E,EMAP,EF,EI,cost,p);
    C.row(e) = p;
    Qit[e] = Q.insert(std::pair<double,int>(cost,e)).first;
  }
  int prev_e = -1;
  bool clean_finish = false;

  while(true)
  {
    if(Q.empty())
    {
      break;
//.........这里部分代码省略.........