C++ eigen::MatrixXi类代码示例

IGL_INLINE void igl::polyvector_field_singularities_from_matchings(
                                                                   const Eigen::PlainObjectBase<DerivedV> &V,
                                                                   const Eigen::PlainObjectBase<DerivedF> &F,
                                                                   const std::vector<bool> &V_border,
                                                                   const std::vector<std::vector<VFType> > &VF,
                                                                   const Eigen::MatrixXi &TT,
                                                                   const Eigen::MatrixXi &E2F,
                                                                   const Eigen::MatrixXi &F2E,
                                                                   const Eigen::PlainObjectBase<DerivedM> &match_ab,
                                                                   const Eigen::PlainObjectBase<DerivedM> &match_ba,
                                                                   Eigen::PlainObjectBase<DerivedS> &singularities,
                                                                   Eigen::PlainObjectBase<DerivedS> &singularity_indices)
{
  igl::polyvector_field_singularities_from_matchings(V, F, V_border, VF, TT, E2F, F2E, match_ab, match_ba, singularities);

  singularity_indices.setZero(singularities.size(), 1);

  //get index from first vector only
  int vector_to_match = 0;
  for (int i =0; i<singularities.size(); ++i)
  {
    int vi = singularities[i];

    // Eigen::VectorXi mvi,fi;
    // igl::polyvector_field_one_ring_matchings(V, F, VF, E2F, F2E, TT, match_ab, match_ba, vi, vector_to_match, mvi, fi);
    Eigen::VectorXi fi;
    Eigen::MatrixXi mvi;
    igl::polyvector_field_one_ring_matchings(V, F, VF, E2F, F2E, TT, match_ab, match_ba, vi, mvi, fi);

    singularity_indices[i] = (mvi(mvi.rows()-1,vector_to_match) - vector_to_match);
  }

}

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

// Converts a representative vector per face in the full set of vectors that describe
// an N-RoSy field
void representative_to_nrosy(
    const Eigen::MatrixXd& V,
    const Eigen::MatrixXi& F,
    const Eigen::MatrixXd& R,
    const int N,
    Eigen::MatrixXd& Y)
{
    using namespace Eigen;
    using namespace std;
    MatrixXd B1, B2, B3;

    igl::local_basis(V,F,B1,B2,B3);

    Y.resize(F.rows()*N,3);
    for (unsigned i=0; i<F.rows(); ++i)
    {
        double x = R.row(i) * B1.row(i).transpose();
        double y = R.row(i) * B2.row(i).transpose();
        double angle = atan2(y,x);

        for (unsigned j=0; j<N; ++j)
        {
            double anglej = angle + 2*M_PI*double(j)/double(N);
            double xj = cos(anglej);
            double yj = sin(anglej);
            Y.row(i*N+j) = xj * B1.row(i) + yj * B2.row(i);
        }
    }
}

void HomogeneousPoint3fIntegralImage::compute(const Eigen::MatrixXi &indices, const HomogeneousPoint3fVector &points) {
  if (cols() != indices.cols() || rows() != indices.rows())
    resize(indices.rows(), indices.cols());
  clear();
  
  HomogeneousPoint3fAccumulator *acc = data();
  const int *pointIndex = indices.data();
  int s = rows() * cols();
  // fill the accumulators with the points
  for (int i=0; i<s; i++, acc++, pointIndex++){
    if (*pointIndex<0)
      continue;
    const HomogeneousPoint3f& point = points[*pointIndex];
    acc->operator += (point);
  }

  // fill by column
  #pragma omp parallel for
  for (int c=0; c<cols(); c++){
    for (int r=1; r<rows(); r++){
      coeffRef(r,c) += coeffRef(r-1,c);
    }
  }

  // fill by row
  #pragma omp parallel for
  for (int r=0; r<rows(); r++){
    for (int c=1; c<cols(); c++){
      coeffRef(r,c) += coeffRef(r,c-1);
    }
  }
}

// Receive a matrix from MATLAB
IGL_INLINE void igl::mlgetmatrix(Engine** mlengine, std::string name, Eigen::MatrixXi& M)
{
  if (*mlengine == 0)
    mlinit(mlengine);
  
  unsigned long m = 0;
  unsigned long n = 0;
  std::vector<double> t;
  
  mxArray *ary = engGetVariable(*mlengine, name.c_str());
  if (ary == NULL)
  {
    m = 0;
    n = 0;
    M = Eigen::MatrixXi(0,0);
  }
  else
  {
    m = mxGetM(ary);
    n = mxGetN(ary);
    M = Eigen::MatrixXi(m,n);
    
    double *pM = mxGetPr(ary);
    
    int c = 0;
    for(int j=0; j<M.cols();++j)
      for(int i=0; i<M.rows();++i)
        M(i,j) = int(pM[c++])-1;
  }
  
  mxDestroyArray(ary);
}

void color_intersections(
  const Eigen::MatrixXd & V,
  const Eigen::MatrixXi & F,
  const Eigen::MatrixXd & U,
  const Eigen::MatrixXi & G,
  Eigen::MatrixXd & C,
  Eigen::MatrixXd & D)
{
  using namespace igl;
  using namespace igl::cgal;
  using namespace Eigen;
  MatrixXi IF;
  const bool first_only = false;
  intersect_other(V,F,U,G,first_only,IF);
  C.resize(F.rows(),3);
  C.col(0).setConstant(0.4);
  C.col(1).setConstant(0.8);
  C.col(2).setConstant(0.3);
  D.resize(G.rows(),3);
  D.col(0).setConstant(0.4);
  D.col(1).setConstant(0.3);
  D.col(2).setConstant(0.8);
  for(int f = 0;f<IF.rows();f++)
  {
    C.row(IF(f,0)) = RowVector3d(1,0.4,0.4);
    D.row(IF(f,1)) = RowVector3d(0.8,0.7,0.3);
  }
}

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

}

IGL_INLINE bool igl::decimate(
  const Eigen::MatrixXd & V,
  const Eigen::MatrixXi & F,
  const size_t max_m,
  Eigen::MatrixXd & U,
  Eigen::MatrixXi & G,
  Eigen::VectorXi & J,
  Eigen::VectorXi & I)
{
  // Original number of faces
  const int orig_m = F.rows();
  // Tracking number of faces
  int m = F.rows();
  typedef Eigen::MatrixXd DerivedV;
  typedef Eigen::MatrixXi DerivedF;
  DerivedV VO;
  DerivedF FO;
  igl::connect_boundary_to_infinity(V,F,VO,FO);
  bool ret = decimate(
    VO,
    FO,
    shortest_edge_and_midpoint,
    max_faces_stopping_condition(m,orig_m,max_m),
    U,
    G,
    J,
    I);
  const Eigen::Array<bool,Eigen::Dynamic,1> keep = (J.array()<orig_m);
  igl::slice_mask(Eigen::MatrixXi(G),keep,1,G);
  igl::slice_mask(Eigen::VectorXi(J),keep,1,J);
  Eigen::VectorXi _1,I2;
  igl::remove_unreferenced(Eigen::MatrixXd(U),Eigen::MatrixXi(G),U,G,_1,I2);
  igl::slice(Eigen::VectorXi(I),I2,1,I);
  return ret;
}

TEST(readOBJ, simple) {
    Eigen::MatrixXd V;
    Eigen::MatrixXi F;
    test_common::load_mesh("cube.obj", V, F);
    ASSERT_EQ(8, V.rows());
    ASSERT_EQ(12, F.rows());
}

IGL_INLINE void igl::covariance_scatter_matrix(
  const Eigen::MatrixXd & V, 
  const Eigen::MatrixXi & F,
  const ARAPEnergyType energy,
  Eigen::SparseMatrix<double>& CSM)
{
  using namespace igl;
  using namespace Eigen;
  // number of mesh vertices
  int n = V.rows();
  assert(n > F.maxCoeff());
  // dimension of mesh
  int dim = V.cols();
  // Number of mesh elements
  int m = F.rows();

  // number of rotations
  int nr;
  switch(energy)
  {
    case ARAP_ENERGY_TYPE_SPOKES:
      nr = n;
      break;
    case ARAP_ENERGY_TYPE_SPOKES_AND_RIMS:
      nr = n;
      break;
    case ARAP_ENERGY_TYPE_ELEMENTS:
      nr = m;
      break;
    default:
      fprintf(
        stderr,
        "covariance_scatter_matrix.h: Error: Unsupported arap energy %d\n",
        energy);
      return;
  }

  SparseMatrix<double> KX,KY,KZ;
  arap_linear_block(V,F,0,energy,KX);
  arap_linear_block(V,F,1,energy,KY);
  SparseMatrix<double> Z(n,nr);
  if(dim == 2)
  {
    CSM = cat(1,cat(2,KX,Z),cat(2,Z,KY)).transpose();
  }else if(dim == 3)
  {
    arap_linear_block(V,F,2,energy,KZ);
    SparseMatrix<double>ZZ(n,nr*2);
    CSM = 
      cat(1,cat(1,cat(2,KX,ZZ),cat(2,cat(2,Z,KY),Z)),cat(2,ZZ,KZ)).transpose();
  }else
  {
    fprintf(
     stderr,
     "covariance_scatter_matrix.h: Error: Unsupported dimension %d\n",
     dim);
    return;
  }

}

IGL_INLINE bool igl::mesh_to_polyhedron(
  const Eigen::MatrixXd & V,
  const Eigen::MatrixXi & F,
  Polyhedron & poly)
{
  typedef typename Polyhedron::HalfedgeDS HalfedgeDS;
  // Postcondition: hds is a valid polyhedral surface.
  CGAL::Polyhedron_incremental_builder_3<HalfedgeDS> B(poly.hds());
  B.begin_surface(V.rows(),F.rows());
  typedef typename HalfedgeDS::Vertex   Vertex;
  typedef typename Vertex::Point Point;
  assert(V.cols() == 3 && "V must be #V by 3");
  for(int v = 0;v<V.rows();v++)
  {
    B.add_vertex(Point(V(v,0),V(v,1),V(v,2)));
  }
  assert(F.cols() == 3 && "F must be #F by 3");
  for(int f=0;f<F.rows();f++)
  {
    B.begin_facet();
    for(int c = 0;c<3;c++)
    {
      B.add_vertex_to_facet(F(f,c));
    }
    B.end_facet();
  }
  if(B.error())
  {
    B.rollback();
    return false;
  }
  B.end_surface();
  return poly.is_valid();
}

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

std::array<Eigen::MatrixXi, 4> FillGroup(const Eigen::MatrixXi& m)
{
  std::array<Eigen::MatrixXi, 4> patterns;
  patterns[0] = m;

  // Found in this awesome answer http://stackoverflow.com/a/3488737/505049
  patterns[1] = m.transpose().colwise().reverse().eval();  // Rotate 90 CW
  patterns[2] = m.transpose().colwise().reverse().transpose().colwise().reverse().eval();  // Rotate 180. Not in the S.O. post
  patterns[3] = m.transpose().rowwise().reverse().eval();  // Rotate 270 CW
  return patterns;
}