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

void update_visualization(igl::opengl::glfw::Viewer & viewer)
{
  using namespace Eigen;
  using namespace std;
  Eigen::Vector4d plane(
    0,0,1,-((1-slice_z)*V.col(2).minCoeff()+slice_z*V.col(2).maxCoeff()));
  MatrixXd V_vis;
  MatrixXi F_vis;
  VectorXi J;
  {
    SparseMatrix<double> bary;
    // Value of plane's implicit function at all vertices
    const VectorXd IV = 
      (V.col(0)*plane(0) + 
        V.col(1)*plane(1) + 
        V.col(2)*plane(2)).array()
      + plane(3);
    igl::marching_tets(V,T,IV,V_vis,F_vis,J,bary);
  }
  VectorXd W_vis;
  igl::slice(W,J,W_vis);
  MatrixXd C_vis;
  // color without normalizing
  igl::parula(W_vis,false,C_vis);


  const auto & append_mesh = [&C_vis,&F_vis,&V_vis](
    const Eigen::MatrixXd & V,
    const Eigen::MatrixXi & F,
    const RowVector3d & color)
  {
    F_vis.conservativeResize(F_vis.rows()+F.rows(),3);
    F_vis.bottomRows(F.rows()) = F.array()+V_vis.rows();
    V_vis.conservativeResize(V_vis.rows()+V.rows(),3);
    V_vis.bottomRows(V.rows()) = V;
    C_vis.conservativeResize(C_vis.rows()+F.rows(),3);
    C_vis.bottomRows(F.rows()).rowwise() = color;
  };
  switch(overlay)
  {
    case OVERLAY_INPUT:
      append_mesh(V,F,RowVector3d(1.,0.894,0.227));
      break;
    case OVERLAY_OUTPUT:
      append_mesh(V,G,RowVector3d(0.8,0.8,0.8));
      break;
    default:
      break;
  }
  viewer.data().clear();
  viewer.data().set_mesh(V_vis,F_vis);
  viewer.data().set_colors(C_vis);
  viewer.data().set_face_based(true);
}

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 update_visualization(igl::viewer::Viewer & viewer)
{
  using namespace Eigen;
  using namespace std;
  Eigen::Vector4d plane(
    0,0,1,-((1-slice_z)*V.col(2).minCoeff()+slice_z*V.col(2).maxCoeff()));
  MatrixXd V_vis;
  MatrixXi F_vis;
  // Extract triangle mesh slice through volume mesh and subdivide nasty
  // triangles
  {
    VectorXi J;
    SparseMatrix<double> bary;
    igl::slice_tets(V,T,plane,V_vis,F_vis,J,bary);
    while(true)
    {
      MatrixXd l;
      igl::edge_lengths(V_vis,F_vis,l);
      l /= (V_vis.colwise().maxCoeff() - V_vis.colwise().minCoeff()).norm();
      const double max_l = 0.03;
      if(l.maxCoeff()<max_l)
      {
        break;
      }
      Array<bool,Dynamic,1> bad = l.array().rowwise().maxCoeff() > max_l;
      MatrixXi F_vis_bad, F_vis_good;
      igl::slice_mask(F_vis,bad,1,F_vis_bad);
      igl::slice_mask(F_vis,(bad!=true).eval(),1,F_vis_good);
      igl::upsample(V_vis,F_vis_bad);
      F_vis = igl::cat(1,F_vis_bad,F_vis_good);
    }
  }

  // Compute signed distance
  VectorXd S_vis;
  {
    VectorXi I;
    MatrixXd N,C;
    // Bunny is a watertight mesh so use pseudonormal for signing
    signed_distance_pseudonormal(V_vis,V,F,tree,FN,VN,EN,EMAP,S_vis,I,C,N);
  }
  // push to [0,1] range
  S_vis.array() = 0.5*(S_vis.array()/max_distance)+0.5;
  MatrixXd C_vis;
  // color without normalizing
  igl::parula(S_vis,false,C_vis);


  const auto & append_mesh = [&C_vis,&F_vis,&V_vis](
    const Eigen::MatrixXd & V,
    const Eigen::MatrixXi & F,
    const RowVector3d & color)
  {
    F_vis.conservativeResize(F_vis.rows()+F.rows(),3);
    F_vis.bottomRows(F.rows()) = F.array()+V_vis.rows();
    V_vis.conservativeResize(V_vis.rows()+V.rows(),3);
    V_vis.bottomRows(V.rows()) = V;
    C_vis.conservativeResize(C_vis.rows()+V.rows(),3);
    C_vis.bottomRows(V.rows()).rowwise() = color;
  };
  if(overlay)
  {
    append_mesh(V,F,RowVector3d(0.8,0.8,0.8));
  }
  viewer.data.clear();
  viewer.data.set_mesh(V_vis,F_vis);
  viewer.data.set_colors(C_vis);
  viewer.core.lighting_factor = overlay;
}

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