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

IGL_INLINE void igl::adjacency_matrix(
  const Eigen::MatrixXi & F, 
  Eigen::SparseMatrix<T>& A)
{
  using namespace std;
  using namespace Eigen;

  typedef Triplet<T> IJV;
  vector<IJV > ijv;
  ijv.reserve(F.size()*2);
  // Loop over faces
  for(int i = 0;i<F.rows();i++)
  {
    // Loop over this face
    for(int j = 0;j<F.cols();j++)
    {
      // Get indices of edge: s --> d
      int s = F(i,j);
      int d = F(i,(j+1)%F.cols());
      ijv.push_back(IJV(s,d,1));
      ijv.push_back(IJV(d,s,1));
    }
  }

  const int n = F.maxCoeff()+1;
  A.resize(n,n);
  switch(F.cols())
  {
    case 3:
      A.reserve(6*(F.maxCoeff()+1));
      break;
    case 4:
      A.reserve(26*(F.maxCoeff()+1));
      break;
  }
  A.setFromTriplets(ijv.begin(),ijv.end());

  // Force all non-zeros to be one

  // Iterate over outside
  for(int k=0; k<A.outerSize(); ++k)
  {
    // Iterate over inside
    for(typename Eigen::SparseMatrix<T>::InnerIterator it (A,k); it; ++it)
    {
      assert(it.value() != 0);
      A.coeffRef(it.row(),it.col()) = 1;
    }
  }
}

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::tetgen::mesh_with_skeleton(
  const Eigen::MatrixXd & V,
  const Eigen::MatrixXi & F,
  const Eigen::MatrixXd & C,
  const Eigen::VectorXi & /*P*/,
  const Eigen::MatrixXi & BE,
  const Eigen::MatrixXi & CE,
  const int samples_per_bone,
  const std::string & tetgen_flags,
  Eigen::MatrixXd & VV,
  Eigen::MatrixXi & TT,
  Eigen::MatrixXi & FF)
{
  using namespace Eigen;
  using namespace igl;
  using namespace std;
  const string eff_tetgen_flags = 
    (tetgen_flags.length() == 0?DEFAULT_TETGEN_FLAGS:tetgen_flags);
  // Collect all edges that need samples:
  MatrixXi BECE = cat(1,BE,CE);
  MatrixXd S;
  // Sample each edge with 10 samples. (Choice of 10 doesn't seem to matter so
  // much, but could under some circumstances)
  sample_edges(C,BECE,samples_per_bone,S);
  // Vertices we'll constrain tet mesh to meet
  MatrixXd VS = cat(1,V,S);
  // Use tetgen to mesh the interior of surface, this assumes surface:
  //   * has no holes
  //   * has no non-manifold edges or vertices
  //   * has consistent orientation
  //   * has no self-intersections
  //   * has no 0-volume pieces
  //writeOBJ("mesh_with_skeleton.obj",VS,F);
  cerr<<"tetgen begin()"<<endl;
  int status = tetrahedralize( VS,F,eff_tetgen_flags,VV,TT,FF);
  cerr<<"tetgen end()"<<endl;
  if(FF.rows() != F.rows())
  {
    // Issue a warning if the surface has changed
    cerr<<"mesh_with_skeleton: Warning: boundary faces != input faces"<<endl;
  }
  if(status != 0)
  {
    cerr<<
      "***************************************************************"<<endl<<
      "***************************************************************"<<endl<<
      "***************************************************************"<<endl<<
      "***************************************************************"<<endl<<
      "* mesh_with_skeleton: tetgen failed. Just meshing convex hull *"<<endl<<
      "***************************************************************"<<endl<<
      "***************************************************************"<<endl<<
      "***************************************************************"<<endl<<
      "***************************************************************"<<endl;
    // If meshing convex hull then use more regular mesh
    status = tetrahedralize(VS,F,"q1.414",VV,TT,FF);
    // I suppose this will fail if the skeleton is outside the mesh
    assert(FF.maxCoeff() < VV.rows());
    if(status != 0)
    {
      cerr<<"mesh_with_skeleton: tetgen failed again."<<endl;
      return false;
    }
  }

  return true;
}

IGL_INLINE void igl::draw_mesh(
  const Eigen::MatrixXd & V,
  const Eigen::MatrixXi & F,
  const Eigen::MatrixXd & N,
  const Eigen::MatrixXi & NF,
  const Eigen::MatrixXd & C,
  const Eigen::MatrixXd & TC,
  const Eigen::MatrixXi & TF,
  const Eigen::MatrixXd & W,
  const GLuint W_index,
  const Eigen::MatrixXi & WI,
  const GLuint WI_index)
{
  using namespace std;
  using namespace Eigen;
  const int rF = F.rows();
  const int cF = F.cols();
  const int cC = C.cols();
  const int rC = C.rows();
  const int cW = W.cols();
  const int rW = W.rows();
  const int rV = V.rows();
  const int rTC = TC.rows();
  const int rTF = TF.rows();
  const int rNF = NF.rows();
  const int rN = N.rows();

  if(F.size() > 0)
  {
    assert(F.maxCoeff() < V.rows());
    assert(V.cols() == 3);
    assert(rC == rV || rC == rF || rC == rF*3 || rC==1 || C.size() == 0);
    assert(C.cols() == 3 || C.size() == 0);
    assert(N.cols() == 3 || N.size() == 0);
    assert(TC.cols() == 2 || TC.size() == 0);
    assert(cF == 3 || cF == 4);
    assert(TF.size() == 0 || TF.cols() == F.cols());
    assert(NF.size() == 0 || NF.cols() == NF.cols());
  }
  if(W.size()>0)
  {
    assert(W.rows() == V.rows());
    assert(WI.rows() == V.rows());
    assert(W.cols() == WI.cols());
  }

  switch(F.cols())
  {
    default:
    case 3:
      glBegin(GL_TRIANGLES);
      break;
    case 4:
      glBegin(GL_QUADS);
      break;
  }
  // loop over faces
  for(int i = 0; i<rF;i++)
  {
    // loop over corners of triangle
    for(int j = 0;j<cF;j++)
    {

      int tc = -1;
      if(rTF != 0)
      {
        tc = TF(i,j);
      } else if(rTC == 1)
      {
        tc = 0;
      }else if(rTC == rV)
      {
        tc = F(i,j);
      }else if(rTC == rF*cF)
      {
        tc = i*cF + j;
      }else if(rTC == rF)
      {
        tc = i;
      }else
      {
        assert(TC.size() == 0);
      }

      // RGB(A)
      Matrix<MatrixXd::Scalar,1,Dynamic> color;
      if(rC == 1)
      {
        color = C.row(0);
      }else if(rC == rV)
      {
        color = C.row(F(i,j));
      }else if(rC == rF*cF)
      {
        color = C.row(i*cF+j);
      }else if(rC == rF)
      {
        color = C.row(i);
      }else
      {
//.........这里部分代码省略.........

bool BF3PointCircle::getRobustCircle(const cvb::CvContourChainCode& contour, const unsigned int maxVotes, const unsigned int maxAccu, const int maxInvalidVotesInSeries, BFCircle& circle) {

	cvb::CvChainCodes::const_iterator it 		= contour.chainCode.begin();
	cvb::CvChainCodes::const_iterator it_beg	= contour.chainCode.begin();
	cvb::CvChainCodes::const_iterator it_end 	= contour.chainCode.end();

	unsigned int x = contour.startingPoint.x;
	unsigned int y = contour.startingPoint.y;

	BFContour bfContour;
	while(it != it_end) {

		bfContour.add(BFCoordinate<int>(static_cast<int>(x),static_cast<int>(y)));

		x += cvb::cvChainCodeMoves[*it][0];
		y += cvb::cvChainCodeMoves[*it][1];

		it++;
	}
	const unsigned int nContourPoints = bfContour.getPixelCount();

	BFRectangle rect = bfContour.getBounds();
	int nRows = bfRound(rect.getHeight());
	int nCols = bfRound(rect.getWidth());
	int x0 = bfRound(rect.getX0());
	int y0 = bfRound(rect.getY0());
	BFCoordinate<int> topLeft(x0,y0);

	// generate 2d histogram for circle center estimation
	Eigen::MatrixXi H = Eigen::MatrixXi::Zero(nRows, nCols);

	unsigned int votes = 0;
	int invalidVotesInSeries = 0;
	while(votes < maxVotes) {

		unsigned int randIndex1 = (rand() % nContourPoints);
		unsigned int randIndex2 = (rand() % nContourPoints);
		while(randIndex2 == randIndex1)
			randIndex2 = (rand() % nContourPoints);
		unsigned int randIndex3 = (rand() % nContourPoints);
		while(randIndex3 == randIndex2 || randIndex3 == randIndex1)
			randIndex3 = (rand() % nContourPoints);

		BFCoordinate<int> c1 = bfContour.getCoordinate(randIndex1) - topLeft;
		BFCoordinate<int> c2 = bfContour.getCoordinate(randIndex2) - topLeft;
		BFCoordinate<int> c3 = bfContour.getCoordinate(randIndex3) - topLeft;

		BFCoordinate<double> center;
		bool validCenter = getCenter(c1,c2,c3,center);

		if(!validCenter) {
			votes--;
			invalidVotesInSeries++;

			if(invalidVotesInSeries > maxInvalidVotesInSeries) {
				return false;
			}
			continue;
		}
		invalidVotesInSeries = 0;

		double cxD = center.getX();
		double cyD = center.getY();

		int cx = bfRound(cxD);
		int cy = bfRound(cyD);

		if(cx < 0 || cy < 0 || cx >= nRows || cy >= nCols) {
			continue;
		}
		else {
			H(cx,cy) += 1;

			if(H(cx,cy) >= static_cast<int>(maxAccu)) {
				break;
			}
		}

		votes++;
	}

	int finalX = 0;
	int finalY = 0;
	H.maxCoeff(&finalX,&finalY);
	finalX += bfRound(x0);
	finalY += bfRound(y0);

	// generate 1d histogram for circle radius estimation
	Eigen::VectorXi K = Eigen::VectorXi::Zero(bfMax(nRows,nCols));

	it = it_beg;
	x = contour.startingPoint.x;
	y = contour.startingPoint.y;

	while(it != it_end) {

		int r = bfRound(sqrt(pow(static_cast<double>(static_cast<int>(x)-finalX),2.0) + pow(static_cast<double>(static_cast<int>(y)-finalY),2.0)));

		if(r < K.rows()) {
			K(r) += 1;
//.........这里部分代码省略.........

void mexFunction(
  int nlhs, mxArray *plhs[], 
  int nrhs, const mxArray *prhs[])
{
  const auto mexErrMsgTxt = [](const bool v ,const char * msg)
  {
    igl::matlab::mexErrMsgTxt(v,msg);
  };
  igl::matlab::MexStream mout;        
  std::streambuf *outbuf = std::cout.rdbuf(&mout);
  mexErrMsgTxt(nrhs >= 4,"Four arguments expected");
  Eigen::MatrixXd V,U0,U,bc;
  Eigen::MatrixXi F;
  Eigen::VectorXi b;
  int iters = 100;
  bool align_guess = true;
  double p = 1e5;

  igl::matlab::parse_rhs_double(prhs+0,V);
  igl::matlab::parse_rhs_index(prhs+1,F);
  igl::matlab::parse_rhs_index(prhs+2,b);
  igl::matlab::parse_rhs_double(prhs+3,bc);

  {
    int i = 4;
    while(i<nrhs)
    {
      mexErrMsgTxt(mxIsChar(prhs[i]),"Parameter names should be strings");
      // Cast to char
      const char * name = mxArrayToString(prhs[i]);
      if(strcmp("P",name) == 0)
      {
        igl::matlab::validate_arg_scalar(i,nrhs,prhs,name);
        igl::matlab::validate_arg_double(i,nrhs,prhs,name);
        p = (double)*mxGetPr(prhs[++i]);
      }else if(strcmp("AlignGuess",name) == 0)
      {
        igl::matlab::validate_arg_scalar(i,nrhs,prhs,name);
        igl::matlab::validate_arg_logical(i,nrhs,prhs,name);
        align_guess = (bool)*mxGetLogicals(prhs[++i]);
      }else if(strcmp("Iters",name) == 0)
      {
        igl::matlab::validate_arg_scalar(i,nrhs,prhs,name);
        igl::matlab::validate_arg_double(i,nrhs,prhs,name);
        iters = (double)*mxGetPr(prhs[++i]);
      }else
      {
        mexErrMsgTxt(false,C_STR("Unknown parameter: "<<name));
      }
      i++;
    }
  }


  {
    Eigen::MatrixXi C;
    igl::vertex_components(F, C);
    mexErrMsgTxt(
      C.maxCoeff() == 0,
      "(V,F) should have exactly 1 connected component");
    const int ec = igl::euler_characteristic(F);
    mexErrMsgTxt(
      ec == 1,
      C_STR("(V,F) should have disk topology (euler characteristic = "<<ec));
    mexErrMsgTxt(
      igl::is_edge_manifold(F),
      "(V,F) should be edge-manifold");
    Eigen::VectorXd A;
    igl::doublearea(V,F,A);
    mexErrMsgTxt(
      (A.array().abs() > igl::EPS<double>()).all(),
      "(V,F) should have non-zero face areas");
  }


  // Initialize with Tutte embedding to disk
  Eigen::VectorXi bnd; 
  Eigen::MatrixXd bnd_uv;
  igl::boundary_loop(F,bnd);
  igl::map_vertices_to_circle(V,bnd,bnd_uv);
  igl::harmonic(V,F,bnd,bnd_uv,1,U0);
  if (igl::flipped_triangles(U0,F).size() != 0) 
  {
    // use uniform laplacian
    igl::harmonic(F,bnd,bnd_uv,1,U0);
  }

  if(align_guess)
  {
    Eigen::MatrixXd X,X0,Y,Y0;
    igl::slice(U0,b,1,X);
    Y = bc;
    Eigen::RowVectorXd Xmean = X.colwise().mean();
    Eigen::RowVectorXd Ymean = Y.colwise().mean();
    Y0 = Y.rowwise()-Ymean;
    X0 = X.rowwise()-Xmean;
    Eigen::MatrixXd I = Eigen::MatrixXd::Identity(X0.cols(),X0.cols())*1e-5;
    Eigen::MatrixXd T = (X0.transpose()*X0+I).inverse()*(X0.transpose()*Y0);
    U0.rowwise() -= Xmean;
    U0 = (U0*T).eval();
//.........这里部分代码省略.........