C++ MatrixXd::colwise方法代码示例

IGL_INLINE void igl::viewer::ViewerCore::get_scale_and_shift_to_fit_mesh(
                                                                 const Eigen::MatrixXd& V,
                                                                 float& zoom,
                                                                 Eigen::Vector3f& shift)
{
  if (V.rows() == 0)
    return;

  Eigen::RowVector3d min_point;
  Eigen::RowVector3d max_point;
  Eigen::RowVector3d centroid;

  if (V.cols() == 3)
  {
    min_point = V.colwise().minCoeff();
    max_point = V.colwise().maxCoeff();
  }
  else if (V.cols() == 2)
  {
    min_point << V.colwise().minCoeff(),0;
    max_point << V.colwise().maxCoeff(),0;
  }
  else
    return;

  centroid = 0.5 * (min_point + max_point);
  shift = -centroid.cast<float>();
  double x_scale = fabs(max_point[0] - min_point[0]);
  double y_scale = fabs(max_point[1] - min_point[1]);
  double z_scale = fabs(max_point[2] - min_point[2]);
  zoom = 2.0 / std::max(z_scale,std::max(x_scale,y_scale));
}

int main(int argc, char * argv[])
{
  using namespace Eigen;
  using namespace std;
  using namespace igl;
  VectorXd D;
  if(!read_triangle_mesh("../shared/beetle.off",V,F))
  {
    cout<<"failed to load mesh"<<endl;
  }
  twod = V.col(2).minCoeff()==V.col(2).maxCoeff();
  bbd = (V.colwise().maxCoeff()-V.colwise().minCoeff()).norm();
  SparseMatrix<double> L,M;
  cotmatrix(V,F,L);
  L = (-L).eval();
  massmatrix(V,F,MASSMATRIX_TYPE_DEFAULT,M);
  const size_t k = 5;
  if(!eigs(L,M,k+1,EIGS_TYPE_SM,U,D))
  {
    cout<<"failed."<<endl;
  }
  U = ((U.array()-U.minCoeff())/(U.maxCoeff()-U.minCoeff())).eval();

  igl::viewer::Viewer viewer;
  viewer.callback_key_down = [&](igl::viewer::Viewer & viewer,unsigned char key,int)->bool
  {
    switch(key)
    {
      default: 
        return false;
      case ' ':
      {
        U = U.rightCols(k).eval();
        // Rescale eigen vectors for visualization
        VectorXd Z = 
          bbd*0.5*U.col(c);
        Eigen::MatrixXd C;
        igl::parula(U.col(c).eval(),false,C);
        c = (c+1)%U.cols();
        if(twod)
        {
          V.col(2) = Z;
        }
        viewer.data.set_mesh(V,F);
        viewer.data.compute_normals();
        viewer.data.set_colors(C);
        return true;
      }
    }
  };
  viewer.callback_key_down(viewer,' ',0);
  viewer.core.show_lines = false;
  viewer.launch();
}

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

  // init mesh
  string filename = "/usr/local/igl/libigl/examples/shared/decimated-knight.obj" ;
  if(argc > 1)
  {
    filename = argv[1];
  }
  if(!read_triangle_mesh(filename,V,F))
  {
    return 1;
  }

  // Compute normals, centroid, colors, bounding box diagonal
  per_face_normals(V,F,N);
  normalize_row_lengths(N,N);
  mid = 0.5*(V.colwise().maxCoeff() + V.colwise().minCoeff());
  bbd =
    (V.colwise().maxCoeff() -
    V.colwise().minCoeff()).maxCoeff();

  // Init glut
  glutInit(&argc,argv);

  if( !TwInit(TW_OPENGL, NULL) )
  {
    // A fatal error occured
    fprintf(stderr, "AntTweakBar initialization failed: %s\n", TwGetLastError());
    return 1;
  }
  // Create a tweak bar
  rebar.TwNewBar("TweakBar");
  rebar.TwAddVarRW("scene_rot", TW_TYPE_QUAT4F, &scene_rot, "");
  rebar.load(REBAR_NAME);

  glutInitDisplayString( "rgba depth double samples>=8 ");
  glutInitWindowSize(glutGet(GLUT_SCREEN_WIDTH)/2.0,glutGet(GLUT_SCREEN_HEIGHT));
  glutCreateWindow("sorted primitives transparency");
  glutDisplayFunc(display);
  glutReshapeFunc(reshape);
  glutKeyboardFunc(key);
  glutMouseFunc(mouse);
  glutMotionFunc(mouse_drag);

  glutMainLoop();

  return 0;
}

IGL_INLINE void igl::viewer::ViewerCore::align_camera_center(
  const Eigen::MatrixXd& V)
{
  if(V.rows() == 0)
    return;

  get_scale_and_shift_to_fit_mesh(V,model_zoom,model_translation);
  // Rather than crash on empty mesh...
  if(V.size() > 0)
  {
    object_scale = (V.colwise().maxCoeff() - V.colwise().minCoeff()).norm();
  }
}

Eigen::MatrixXd localWeighting( const Eigen::MatrixXd &W, bool isFull, bool isSymmetric)
{
	int n = W.rows();
	double Ls = (W.count()-n)/2; //count number of edges of the graph

	Eigen::MatrixXd C = Eigen::MatrixXd::Zero(n,n);
	double Ws = 0.5*W.sum();

	Eigen::VectorXd D = W.colwise().sum();

	if (isFull)
	{
		if (isSymmetric)
		{
			computeC_symmetric_full(const_cast<Eigen::MatrixXd &>(W),C,D,n);
		}
		else
		{
			// this is a trick to ensure vectorizatoin and no cache misses! some tranpositions have to be made though
			const_cast<Eigen::MatrixXd &>(W).transposeInPlace();
			computeC_symmetric_full(const_cast<Eigen::MatrixXd &>(W),C,D,n);
			const_cast<Eigen::MatrixXd &>(W).transposeInPlace();
			C.transposeInPlace();
			// the original code use vertical access to rows, but is slower
			//compute_C_asymmetric_full(const_cast<Eigen::MatrixXd &>(W),C,D,n);
		}
	}
	else
	{
		if (isSymmetric)
		{
			computeC_symmetric_sparse(const_cast<Eigen::MatrixXd &>(W),C,D,n);
		}
		else
		{
			compute_C_asymmetric_sparse(const_cast<Eigen::MatrixXd &>(W),C,D,n);
		}
	}

	Eigen::MatrixXd G = ((Ls/Ws)*W).cwiseProduct(C);

	Eigen::VectorXd DG = G.colwise().sum();

	for (int i=0; i<n; i++)
	{
		G.row(i)/=DG(i);
	}
	
	return G;
}

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

  // init mesh
  string filename = "../shared/decimated-knight.obj";
  if(argc < 2)
  {
    cerr<<"Usage:"<<endl<<"    ./example input.obj"<<endl;
    cout<<endl<<"Opening default mesh..."<<endl;
  }else
  {
    // Read and prepare mesh
    filename = argv[1];
  }
  if(!read_triangle_mesh(filename,V,F))
  {
    return 1;
  }
  // Compute normals, centroid, colors, bounding box diagonal
  per_face_normals(V,F,N);
  normalize_row_lengths(N,N);
  mean = V.colwise().mean();
  C.resize(F.rows(),3);
  init_C();
  bbd =
    (V.colwise().maxCoeff() -
    V.colwise().minCoeff()).maxCoeff();

  // Init embree
  ei.init(V.cast<float>(),F.cast<int>());

  // Init glut
  glutInit(&argc,argv);
  glutInitDisplayString( "rgba depth double samples>=8 ");
  glutInitWindowSize(glutGet(GLUT_SCREEN_WIDTH)/2.0,glutGet(GLUT_SCREEN_HEIGHT));
  glutCreateWindow("embree");
  glutDisplayFunc(display);
  glutReshapeFunc(reshape);
  glutKeyboardFunc(key);
  glutMouseFunc(mouse);
  glutMotionFunc(mouse_drag);
  glutPassiveMotionFunc(mouse_move);
  glutMainLoop();
  return 0;
}

Eigen::MatrixXd caculateAlgebraicDistance(const Eigen::SparseMatrix<int,Eigen::RowMajor>& incidenceMatrix,
                                          const Eigen::SparseMatrix<double,Eigen::RowMajor>& updatedMatrix)
{
  Eigen::MatrixXd inc = Eigen::MatrixXd(incidenceMatrix);
  Eigen::MatrixXd update = Eigen::MatrixXd(updatedMatrix);
  int edgeNumber = inc.cols();
  int verticeNumber = inc.rows();
  Eigen::MatrixXd algebraicDistanceForEdge(edgeNumber,1);
  Eigen::VectorXd sumOfcolsVector(edgeNumber);
  sumOfcolsVector = inc.colwise().sum();
  for(int i = 0; i<edgeNumber;++i){
    int verticeNumberForEdge = sumOfcolsVector[i];
    Eigen::MatrixXd allEdgeMatrix(update.rows(),verticeNumberForEdge);
    int g=0;
    for(int j=0;j<verticeNumber;++j){
        if(inc(j,i)){
          ++g;
          for(int v = 0; v<update.rows();++v){
              allEdgeMatrix(v,g-1) = update(v,j);
          }
        }
    }

    if(g!=0){
      Eigen::MatrixXd maxValueOfEdge = getMaxDiffForeachEdge(allEdgeMatrix);
      maxValueOfEdge = maxValueOfEdge.transpose() * maxValueOfEdge;
      double edgeDistance = maxValueOfEdge.sum();
      edgeDistance = sqrt(edgeDistance);
      algebraicDistanceForEdge(i,0) = edgeDistance;

    }

  }
  return algebraicDistanceForEdge;
}

bool CGEOM::generate_scene_trans
( const SceneGeneratorOptions& sc_opts,
  Eigen::MatrixXd& mP3D,
  Eigen::MatrixXd& mMeasT,
  Eigen::MatrixXd& mMeasN,
  Eigen::Matrix3d& mR,
  Eigen::Vector3d& mt
) {
    if( !generate_scene( sc_opts,
                         mP3D, mMeasT, mMeasN
                       ) ) {
        return false;
    }

    // Create random transform
    mt = mP3D.rowwise().mean();
    const double drotx = rand_range_d( -45., 45. )*3.14159/180.;
    const double droty = rand_range_d( -45., 45. )*3.14159/180.;
    const double drotz = rand_range_d( -45., 45. )*3.14159/180.;
#if 1
    mR =
        ( CEIGEN::skew_rot<Eigen::Matrix3d,double>( drotx, 0., 0. ) +
          CEIGEN::skew_rot<Eigen::Matrix3d,double>( 0., droty , 0.) +
          CEIGEN::skew_rot<Eigen::Matrix3d,double>( 0., 0., drotz ) ).exp();
#else
    mR = Eigen::Matrix3d::Identity();
#endif

    mP3D.colwise() -= mt;
    mP3D = mR.transpose() * mP3D;

    return true;
}

IGL_INLINE void igl::viewer::ViewerCore::get_scale_and_shift_to_fit_mesh(
  const Eigen::MatrixXd& V,
  const Eigen::MatrixXi& F,
  float& zoom,
  Eigen::Vector3f& shift)
{
  if (V.rows() == 0)
    return;

  //Eigen::SparseMatrix<double> M;
  //igl::massmatrix(V,F,igl::MASSMATRIX_TYPE_VORONOI,M);
  //const auto & MV = M*V;
  //Eigen::RowVector3d centroid  = MV.colwise().sum()/M.diagonal().sum();
  Eigen::MatrixXd BC;
  igl::barycenter(V,F,BC);
  Eigen::RowVector3d min_point = BC.colwise().minCoeff();
  Eigen::RowVector3d max_point = BC.colwise().maxCoeff();
  Eigen::RowVector3d centroid  = 0.5*(min_point + max_point);

  shift = -centroid.cast<float>();
  double x_scale = fabs(max_point[0] - min_point[0]);
  double y_scale = fabs(max_point[1] - min_point[1]);
  double z_scale = fabs(max_point[2] - min_point[2]);
  zoom = 2.0 / std::max(z_scale,std::max(x_scale,y_scale));
}

IGL_INLINE void igl::viewer::ViewerCore::get_scale_and_shift_to_fit_mesh(
  const Eigen::MatrixXd& V,
  const Eigen::MatrixXi& F,
  float& zoom,
  Eigen::Vector3f& shift)
{
  if (V.rows() == 0)
    return;

  Eigen::MatrixXd BC;
  if (F.rows() <= 1)
    BC = V;
  else
    igl::barycenter(V,F,BC);

  Eigen::RowVector3d min_point = BC.colwise().minCoeff();
  Eigen::RowVector3d max_point = BC.colwise().maxCoeff();
  Eigen::RowVector3d centroid  = 0.5*(min_point + max_point);

  shift = -centroid.cast<float>();
  double x_scale = fabs(max_point[0] - min_point[0]);
  double y_scale = fabs(max_point[1] - min_point[1]);
  double z_scale = fabs(max_point[2] - min_point[2]);
  zoom = 2.0 / std::max(z_scale,std::max(x_scale,y_scale));
}

int main(int argc, char *argv[])
{
  using namespace Eigen;
  using namespace std;
  igl::readOFF("../shared/decimated-knight.off",V,F);
  U=V;
  igl::readDMAT("../shared/decimated-knight-selection.dmat",S);

  // vertices in selection
  igl::colon<int>(0,V.rows()-1,b);
  b.conservativeResize(stable_partition( b.data(), b.data()+b.size(), 
   [](int i)->bool{return S(i)>=0;})-b.data());
  // Centroid
  mid = 0.5*(V.colwise().maxCoeff() + V.colwise().minCoeff());
  // Precomputation
  arap_data.max_iter = 100;
  igl::arap_precomputation(V,F,V.cols(),b,arap_data);

  // Set color based on selection
  MatrixXd C(F.rows(),3);
  RowVector3d purple(80.0/255.0,64.0/255.0,255.0/255.0);
  RowVector3d gold(255.0/255.0,228.0/255.0,58.0/255.0);
  for(int f = 0;f<F.rows();f++)
  {
    if( S(F(f,0))>=0 && S(F(f,1))>=0 && S(F(f,2))>=0)
    {
      C.row(f) = purple;
    }else
    {
      C.row(f) = gold;
    }
  }

  // Plot the mesh with pseudocolors
  igl::viewer::Viewer viewer;
  viewer.data.set_mesh(U, F);
  viewer.data.set_colors(C);
  viewer.callback_pre_draw = &pre_draw;
  viewer.callback_key_down = &key_down;
  viewer.core.is_animating = false;
  viewer.core.animation_max_fps = 30.;
  cout<<
    "Press [space] to toggle animation"<<endl;
  viewer.launch();
}

std::pair<Eigen::MatrixXd::Index, double> dsearch(Eigen::Vector3d p, Eigen::MatrixXd positions)
{
	Eigen::MatrixXd::Index index;
	// find nearest neighbour
	(positions.colwise() - p).colwise().squaredNorm().minCoeff(&index);
	double d = (positions.col(index) - p).norm();

	return std::make_pair(index , d);
}

void init_relative()
{
  using namespace Eigen;
  using namespace igl;
  per_face_normals(V,F,N);
  Vmax = V.colwise().maxCoeff();
  Vmin = V.colwise().minCoeff();
  Vmid = 0.5*(Vmax + Vmin);
  bbd = (Vmax-Vmin).norm();
}

int FeatureTransformationEstimator::consensus3D(Eigen::MatrixXd P, Eigen::MatrixXd Q, Eigen::Isometry3d T, double thresh, Eigen::Array<bool, 1, Eigen::Dynamic> &consensusSet)
{
    int consensus = 0;

    P = T * P.colwise().homogeneous();
    Eigen::MatrixXd norms = (P - Q).colwise().norm();
    consensusSet = norms.array() < thresh;
    consensus = consensusSet.count();

    return consensus;
}

void init_relative()
{
  using namespace Eigen;
  using namespace igl;
  using namespace std;
  per_face_normals(V,F,N);
  const auto Vmax = V.colwise().maxCoeff();
  const auto Vmin = V.colwise().minCoeff();
  Vmid = 0.5*(Vmax + Vmin);
  centroid(V,F,Vcen);
  bbd = (Vmax-Vmin).norm();
  camera.push_away(2);
}