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

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

  // Load a quad mesh generated by a conjugate field
  igl::readOFF(TUTORIAL_SHARED_PATH "/inspired_mesh_quads_Conjugate.off", VQC, FQC);

  // Convert it in a triangle mesh
  FQCtri.resize(2*FQC.rows(), 3);
  FQCtri <<  FQC.col(0),FQC.col(1),FQC.col(2),
             FQC.col(2),FQC.col(3),FQC.col(0);
  igl::slice( VQC, FQC.col(0).eval(), 1, PQC0);
  igl::slice( VQC, FQC.col(1).eval(), 1, PQC1);
  igl::slice( VQC, FQC.col(2).eval(), 1, PQC2);
  igl::slice( VQC, FQC.col(3).eval(), 1, PQC3);

  // Planarize it
  igl::planarize_quad_mesh(VQC, FQC, 100, 0.005, VQCplan);

  // Convert the planarized mesh to triangles
  igl::slice( VQCplan, FQC.col(0).eval(), 1, PQC0plan);
  igl::slice( VQCplan, FQC.col(1).eval(), 1, PQC1plan);
  igl::slice( VQCplan, FQC.col(2).eval(), 1, PQC2plan);
  igl::slice( VQCplan, FQC.col(3).eval(), 1, PQC3plan);

  // Launch the viewer
  igl::viewer::Viewer viewer;
  key_down(viewer,'2',0);
  viewer.data.invert_normals = true;
  viewer.data.show_lines = false;
  viewer.callback_key_down = &key_down;
  viewer.launch();
}

void PixelMapper::mergeProjections(Eigen::MatrixXf& depthImage, Eigen::MatrixXi& indexImage,
				   Eigen::MatrixXf* depths, Eigen::MatrixXi* indices, int numImages){
  assert (numImages>0);
  int rows=depths[0].rows();
  int cols=depths[0].cols();
  depthImage.resize(indexImage.rows(), indexImage.cols());
  depthImage.fill(std::numeric_limits<float>::max());
  indexImage.resize(rows, cols);
  indexImage.fill(-1);

  #pragma omp parallel for
  for (int c=0; c<cols; c++){
    int* destIndexPtr = &indexImage.coeffRef(0,c);
    float* destDepthPtr = &depthImage.coeffRef(0,c);
    int* srcIndexPtr[numImages];
    float* srcDepthPtr[numImages];
    for (int i=0; i<numImages; i++){
      srcIndexPtr[i] = &indices[i].coeffRef(0,c);
      srcDepthPtr[i] = &depths[i].coeffRef(0,c);
    }
    for (int r=0; r<rows; r++){
      for (int i=0; i<numImages; i++){
	if (*destDepthPtr>*srcDepthPtr[i]){
	  *destDepthPtr = *srcDepthPtr[i];
	  *destIndexPtr = *srcIndexPtr[i];
	}
	srcDepthPtr[i]++;
	srcIndexPtr[i]++;
      }
      destDepthPtr++;
      destIndexPtr++;
    }
  }

}

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

  // Create the boundary of a square
  V.resize(8,2);
  E.resize(8,2);
  H.resize(1,2);

  V << -1,-1, 1,-1, 1,1, -1, 1,
       -2,-2, 2,-2, 2,2, -2, 2;

  E << 0,1, 1,2, 2,3, 3,0,
       4,5, 5,6, 6,7, 7,4;

  H << 0,0;

  // Triangulate the interior
  igl::triangle::triangulate(V,E,H,"a0.005q",V2,F2);

  // Plot the generated mesh
  igl::viewer::Viewer viewer;
  viewer.data.set_mesh(V2,F2);
  viewer.launch();
}

//=============================================================================
int readMatrixXi(Eigen::MatrixXi &_matrix, const char* _filename)
{
	std::ifstream ifs(_filename);

	int rows;
	ifs >> rows;

	int cols;
	ifs >> cols;

	if (rows < 1 || cols < 1)
	{
		std::cout << "Invalid matrix size: (" << rows << ", " << cols << ")." 
			<< std::endl;

		return -1;
	}

	_matrix.resize(rows, cols);

	for (int i = 0; i < rows; i++) {
		for (int j = 0; j < cols; j++)
		{
			ifs >> _matrix(i, j);
		}
	}

	return 0;
}

// Parse right hand side arguments for a matlab mex function.
//
// Inputs:
//   nrhs  number of right hand side arguments
//   prhs  pointer to right hand side arguments
// Outputs:
//   V  n by dim list of mesh vertex positions
//   F  m by dim list of mesh face indices
//   s  1 by dim bone source vertex position
//   d  1 by dim bone dest vertex position
// "Throws" matlab errors if dimensions are not sane.
void parse_rhs(
  const int nrhs, 
  const mxArray *prhs[], 
  Eigen::MatrixXd & V,
  Eigen::MatrixXi & F,
  Eigen::VectorXd & s,
  Eigen::VectorXd & d)
{
  using namespace std;
  if(nrhs < 4)
  {
    mexErrMsgTxt("nrhs < 4");
  }
  const int dim = mxGetN(prhs[0]);
  if(dim != 3)
  {
    mexErrMsgTxt("Mesh vertex list must be #V by 3 list of vertex positions");
  }
  if(dim != (int)mxGetN(prhs[1]))
  {
   mexErrMsgTxt("Mesh facet size must equal dimension");
  }
  if(dim != (int)mxGetN(prhs[2]))
  {
   mexErrMsgTxt("Source dim must equal vertex dimension");
  }
  if(dim != (int)mxGetN(prhs[3]))
  {
   mexErrMsgTxt("Dest dim must equal vertex dimension");
  }
  // set number of mesh vertices
  const int n = mxGetM(prhs[0]);
  // set vertex position pointers
  double * Vp = mxGetPr(prhs[0]);
  // set number of faces
  const int m = mxGetM(prhs[1]);
  // set face index list pointer
  double * Fp = mxGetPr(prhs[1]);
  // set source and dest pointers
  double * sp = mxGetPr(prhs[2]);
  double * dp = mxGetPr(prhs[3]);
  // resize output to transpose
  V.resize(n,dim);
  copy(Vp,Vp+n*dim,V.data());
  // resize output to transpose
  F.resize(m,dim);
  // Q: Is this doing a cast?
  // A: Yes.
  copy(Fp,Fp+m*dim,F.data());
  // http://stackoverflow.com/a/4461466/148668
  transform(F.data(),F.data()+m*dim,F.data(),
    bind2nd(std::plus<double>(),-1.0));
  // resize output to transpose
  s.resize(dim);
  copy(sp,sp+dim,s.data());
  d.resize(dim);
  copy(dp,dp+dim,d.data());
}

IGL_INLINE void igl::exterior_edges(
  const Eigen::MatrixXi & F,
  Eigen::MatrixXi & E)
{
  using namespace Eigen;
  using namespace std;
  assert(F.cols() == 3);
  const size_t m = F.rows();
  MatrixXi all_E,sall_E,sort_order;
  // Sort each edge by index
  all_edges(F,all_E);
  sort(all_E,2,true,sall_E,sort_order);
  // Find unique edges
  MatrixXi uE;
  VectorXi IA,EMAP;
  unique_rows(sall_E,uE,IA,EMAP);
  VectorXi counts = VectorXi::Zero(uE.rows());
  for(size_t a = 0;a<3*m;a++)
  {
    counts(EMAP(a)) += (sort_order(a)==0?1:-1);
  }

  E.resize(all_E.rows(),2);
  {
    int e = 0;
    const size_t nue = uE.rows();
    // Append each unique edge with a non-zero amount of signed occurances
    for(size_t ue = 0; ue<nue; ue++)
    {
      const int count = counts(ue);
      size_t i,j;
      if(count == 0)
      {
        continue;
      }else if(count < 0)
      {
        i = uE(ue,1);
        j = uE(ue,0);
      }else if(count > 0)
      {
        i = uE(ue,0);
        j = uE(ue,1);
      }
      // Append edge for every repeated entry
      const int abs_count = abs(count);
      for(size_t k = 0;k<abs_count;k++)
      {
        E(e,0) = i;
        E(e,1) = j;
        e++;
      }
    }
    E.conservativeResize(e,2);
  }
}

 // construct
 WorldParamsPrior(const std::vector<distrib::MultiGaussian>& ctrlpts_, const std::vector<Eigen::MatrixXi>& ctrlptMasks_,
                  const std::vector<Eigen::MatrixXd>& ctrlptMins_, const std::vector<Eigen::MatrixXd>& ctrlptMaxs_,
                  const Eigen::VectorXd& ctrlptCoupledSds_, const Eigen::VectorXd& ctrlptUncoupledSds_,
                  const std::vector<distrib::MultiGaussian>& properties_, const std::vector<Eigen::VectorXi>& propMasks_,
                  const std::vector<Eigen::VectorXd>& propMins_, const std::vector<Eigen::VectorXd>& propMaxs_,
                  const std::vector<BoundaryClass>& classes_)
     : ctrlptMasks(ctrlptMasks_),
       ctrlptMins(ctrlptMins_),
       ctrlptMaxs(ctrlptMaxs_),
       ctrlptCoupledSds(ctrlptCoupledSds_),
       ctrlptUncoupledSds(ctrlptUncoupledSds_),
       ctrlptPrior(ctrlpts_),
       propMasks(propMasks_),
       propMins(propMins_),
       propMaxs(propMaxs_),
       propertyPrior(properties_),
       classes(classes_)
 {
   for (uint l = 0; l < propMasks.size(); l++)
   {
     for (uint p = 0; p < propMasks_[l].size(); p++)
     {
       if (!propMasks_[l](p))
       {
         propertyPrior[l].sigma.row(p).setZero();
         propertyPrior[l].sigma.col(p).setZero();
         propertyPrior[l].sigma(p, p) = 1;
       }
     }
   }
   for (uint l = 0; l < ctrlptMasks_.size(); l++)
   {
     Eigen::MatrixXi masks = ctrlptMasks_[l];
     masks.resize(masks.rows() * masks.cols(), 1);
     for (uint p = 0; p < masks.rows(); p++)
     {
       if (!masks(p))
       {
         ctrlptPrior[l].sigma.row(p).setZero();
         ctrlptPrior[l].sigma.col(p).setZero();
         ctrlptPrior[l].sigma(p, p) = 1;
       }
     }
   }
   WorldParams minParams; 
   WorldParams maxParams;
   minParams.rockProperties = propMins;
   maxParams.rockProperties = propMaxs;
   minParams.controlPoints = ctrlptMins;
   maxParams.controlPoints = ctrlptMaxs;
   thetaMin = deconstruct(minParams);
   thetaMax = deconstruct(maxParams);
   VLOG(1) << "Theta min:" << thetaMin.transpose();
   VLOG(1) << "Theta max:" << thetaMax.transpose();
 }

void parse_rhs(
  const int nrhs, 
  const mxArray *prhs[], 
  Eigen::MatrixXd & V,
  Eigen::MatrixXi & F,
  Eigen::MatrixXd & P,
  Eigen::MatrixXd & N,
  int & num_samples)
{
  using namespace std;
  if(nrhs < 5)
  {
    mexErrMsgTxt("nrhs < 5");
  }

  const int dim = mxGetN(prhs[0]);
  if(dim != 3)
  {
    mexErrMsgTxt("Mesh vertex list must be #V by 3 list of vertex positions");
  }
  if(dim != (int)mxGetN(prhs[1]))
  {
   mexErrMsgTxt("Mesh facet size must be 3");
  }
  if(mxGetN(prhs[2]) != dim)
  {
    mexErrMsgTxt("Point list must be #P by 3 list of origin locations");
  }
  if(mxGetN(prhs[3]) != dim)
  {
    mexErrMsgTxt("Normal list must be #P by 3 list of origin normals");
  }
  if(mxGetN(prhs[4]) != 1 || mxGetM(prhs[4]) != 1)
  {
    mexErrMsgTxt("Number of samples must be scalar.");
  }


  V.resize(mxGetM(prhs[0]),mxGetN(prhs[0]));
  copy(mxGetPr(prhs[0]),mxGetPr(prhs[0])+V.size(),V.data());
  F.resize(mxGetM(prhs[1]),mxGetN(prhs[1]));
  copy(mxGetPr(prhs[1]),mxGetPr(prhs[1])+F.size(),F.data());
  F.array() -= 1;
  P.resize(mxGetM(prhs[2]),mxGetN(prhs[2]));
  copy(mxGetPr(prhs[2]),mxGetPr(prhs[2])+P.size(),P.data());
  N.resize(mxGetM(prhs[3]),mxGetN(prhs[3]));
  copy(mxGetPr(prhs[3]),mxGetPr(prhs[3])+N.size(),N.data());
  if(*mxGetPr(prhs[4]) != (int)*mxGetPr(prhs[4]))
  {
    mexErrMsgTxt("Number of samples should be non negative integer.");
  }
  num_samples = (int) *mxGetPr(prhs[4]);
}

void PointProjector::projectIntervals(Eigen::MatrixXi &intervalImage, 
				      const Eigen::MatrixXf &depthImage, 
				      const float worldRadius) const{
  intervalImage.resize(depthImage.rows(), depthImage.cols());
  int cpix=0;
  for (int c=0; c<depthImage.cols(); c++){
    const float* f = &depthImage(0,c);
    int* i =&intervalImage(0,c);
    for (int r=0; r<depthImage.rows(); r++, f++, i++){
      *i=projectInterval(r,c,*f, worldRadius);
      cpix++;
    }
  }
}

IGL_INLINE void igl::copyleft::cgal::polyhedron_to_mesh(
  const Polyhedron & poly,
  Eigen::MatrixXd & V,
  Eigen::MatrixXi & F)
{
  using namespace std;
  V.resize(poly.size_of_vertices(),3);
  F.resize(poly.size_of_facets(),3);
  typedef typename Polyhedron::Vertex_const_iterator Vertex_iterator;
  std::map<Vertex_iterator,size_t> vertex_to_index;
  {
    size_t v = 0;
    for(
      typename Polyhedron::Vertex_const_iterator p = poly.vertices_begin();
      p != poly.vertices_end();
      p++)
    {
      V(v,0) = p->point().x();
      V(v,1) = p->point().y();
      V(v,2) = p->point().z();
      vertex_to_index[p] = v;
      v++;
    }
  }
  {
    size_t f = 0;
    for(
      typename Polyhedron::Facet_const_iterator facet = poly.facets_begin();
      facet != poly.facets_end();
      ++facet)
    {
      typename Polyhedron::Halfedge_around_facet_const_circulator he = 
        facet->facet_begin();
      // Facets in polyhedral surfaces are at least triangles.
      assert(CGAL::circulator_size(he) == 3 && "Facets should be triangles");
      size_t c = 0;
      do {
        //// This is stooopidly slow
        // F(f,c) = std::distance(poly.vertices_begin(), he->vertex());
        F(f,c) = vertex_to_index[he->vertex()];
        c++;
      } while ( ++he != facet->facet_begin());
      f++;
    }
  }
}

/**
 * @brief Searches in the graph closest points to origin and target
 * 
 * @param origin ...
 * @param target ...
 * @param originVertex to return selected vertex
 * @param targetVertex ...
 * @return void
 */
void PlannerPRM::searchClosestPoints(const QVec& origin, const QVec& target, Vertex& originVertex, Vertex& targetVertex)
{
	qDebug() << __FUNCTION__ << "Searching from " << origin << "and " << target;

	//prepare the query
	Eigen::MatrixXi indices;
	Eigen::MatrixXf distsTo;
	Eigen::MatrixXf query(3,2);
	indices.resize(1, query.cols());
	distsTo.resize(1, query.cols());
	query(0,0) = origin.x();query(1,0) = origin.y();query(2,0) = origin.z();
	query(0,1) = target.x();query(1,1) = target.y();query(2,1) = target.z();

	nabo->knn(query, indices, distsTo, 1);

	originVertex = vertexMap.value(indices(0,0));
	targetVertex = vertexMap.value(indices(0,1));

 	qDebug() << __FUNCTION__ << "Closest point to origin is at" << data(0,indices(0,0)) << data(1,indices(0,0)) << data(2,indices(0,0)) << " and corresponds to " << graph[originVertex].pose;
 	qDebug() << __FUNCTION__ << "Closest point to target is at" << data(0,indices(0,1)) << data(1,indices(0,1)) << data(2,indices(0,1)) << " and corresponds to " << graph[targetVertex].pose;

}

IGL_INLINE void igl::triangle_fan(
  const Eigen::MatrixXi & E,
  Eigen::MatrixXi & cap)
{
  using namespace std;
  using namespace Eigen;

  // Handle lame base case
  if(E.size() == 0)
  {
    cap.resize(0,E.cols()+1);
    return;
  }
  // "Triangulate" aka "close" the E trivially with facets
  // Note: in 2D we need to know if E endpoints are incoming or
  // outgoing (left or right). Thus this will not work.
  assert(E.cols() == 2);
  // Arbitrary starting vertex
  //int s = E(int(((double)rand() / RAND_MAX)*E.rows()),0);
  int s = E(rand()%E.rows(),0);
  vector<vector<int> >  lcap;
  for(int i = 0;i<E.rows();i++)
  {
    // Skip edges incident on s (they would be zero-area)
    if(E(i,0) == s || E(i,1) == s)
    {
      continue;
    }
    vector<int> e(3);
    e[0] = s;
    e[1] = E(i,0);
    e[2] = E(i,1);
    lcap.push_back(e);
  }
  list_to_matrix(lcap,cap);
}

void PointProjector::unProject(HomogeneousPoint3fVector &points,
			       Eigen::MatrixXi &indexImage,
			       const Eigen::MatrixXf &depthImage) const {
  points.resize(depthImage.rows()*depthImage.cols());
  int count = 0;
  indexImage.resize(depthImage.rows(), depthImage.cols());
  HomogeneousPoint3f* point = &points[0];
  int cpix=0;
  for (int c=0; c<depthImage.cols(); c++){
    const float* f = &depthImage(0,c);
    int* i =&indexImage(0,c);
    for (int r=0; r<depthImage.rows(); r++, f++, i++){
      if (!unProject(*point, r,c,*f)){
	*i=-1;
	continue;
      }
      point++;
      cpix++;
      *i=count;
      count++;
    }
  }
  points.resize(count);
}

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

//.........这里部分代码省略.........
      //
      case1In(dd, T.row(i), inside, outside, identical, inside_t, outside_t);
    } else if (inside.size() == 3 && outside.size() == 1) {
      // (these are colored orange)
      times[3]++;
      // (C) takes care of:
      //   inside: 3 outside: 1
      //
      case3In1Out(dd, T.row(i), inside, outside, identical, inside_t, outside_t);
    } else if (inside.size() == 2 && outside.size() >= 1) {
      // (these are colored red)
      times[4]++;
      // (D) takes care of:
      //   inside: 2 outside: 1 identical: 1
      //   inside: 2 outside: 2 identical: 0
      //
      case2In2Out(dd, T.row(i), inside, outside, identical, inside_t, outside_t);
    } else {
      times[5]++;
      fprintf(stderr, "WARN: marching tets found something weird, with in:%lu out:%lu\n",
              inside.size(), outside.size());
    }
  }

  printf("Finished marching tets with usages:\n");
  for (int i = 0; i < 6; ++i) {
    printf("  %d: %d\n", i, times[i]);
  }
  printf("how_much is %lf and EPS is %lf\n", min_how_much, GLOBAL::EPS);
  printf("     max is %lf\n", max_how_much);
  printf("Num equal is %d\n", numEq);

  // Copy verts
  NV.resize(V.rows() + new_verts.size(), 3);
  for (int i = 0; i < V.rows(); ++i) {
    NV.row(i) = V.row(i);
  }
  for (int i = 0; i < new_verts.size(); ++i) {
    NV.row(i + V.rows()) = new_verts[i];
  }
  // Set I
  I.resize(NV.rows());
  for (int i = 0; i < I.rows(); ++i) {
    if (i < V.rows()) {
      I(i) = i;
    } else {
      I(i) = -1;
    }
  }
  Eigen::VectorXi facesMarkers;
  facesMarkers.resize(faces.size());
  // Copy faces
  NF.resize(faces.size(), 3);
  for (int i = 0; i < faces.size(); ++i) {
    NF.row(i) = faces[i];
    facesMarkers(i) = faces_markers[i];
  }
  
  Eigen::MatrixXd newV;
  Eigen::MatrixXi newF;
  Eigen::VectorXi SVJ, SVI, I2;
  // Helpers::viewTriMesh(NV, NF, facesMarkers);
  //igl::writeOFF("offset_mesh.off", NV, NF)
  Helpers::writeMeshWithMarkers("offset_mesh", NV, NF, facesMarkers);

  /*