You can not select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
122 lines
3.7 KiB
122 lines
3.7 KiB
#include <igl/readOBJ.h>
|
|
#include <igl/cotmatrix.h>
|
|
#include <igl/massmatrix.h>
|
|
#include <igl/edges.h>
|
|
#include <igl/find.h>
|
|
#include <igl/min_quad_with_fixed.h>
|
|
#include <igl/polygons_to_triangles.h>
|
|
#include <igl/opengl/glfw/Viewer.h>
|
|
#include <iostream>
|
|
|
|
int main(int argc, char *argv[])
|
|
{
|
|
using namespace Eigen;
|
|
using namespace std;
|
|
|
|
// Load a mesh in OBJ format
|
|
Eigen::MatrixXd OV,V;
|
|
Eigen::VectorXi I,C;
|
|
igl::readOBJ(
|
|
argc<=1?TUTORIAL_SHARED_PATH "/cylinder.obj" :argv[1],V,I,C);
|
|
// center
|
|
V.rowwise() -= V.colwise().mean();
|
|
OV = V;
|
|
// Convert polygon representation to triangles
|
|
Eigen::MatrixXi F;
|
|
Eigen::VectorXi J;
|
|
igl::polygons_to_triangles(I,C,F,J);
|
|
Eigen::SparseMatrix<double> pL,pM,pP;
|
|
igl::cotmatrix(V,I,C,pL,pM,pP);
|
|
Eigen::SparseMatrix<double> tL,tM;
|
|
igl::cotmatrix(V,F,tL);
|
|
igl::massmatrix(V,F,igl::MASSMATRIX_TYPE_DEFAULT,tM);
|
|
const double bbd = (V.colwise().maxCoeff()- V.colwise().minCoeff()).norm();
|
|
igl::opengl::glfw::Viewer vr;
|
|
vr.data_list[0].set_mesh(V,F);
|
|
vr.append_mesh();
|
|
vr.selected_data_index = 0;
|
|
vr.data_list[0].set_face_based(true);
|
|
Eigen::MatrixXi E;
|
|
igl::edges(I,C,E);
|
|
bool show_edges = true;
|
|
bool use_poly = true;
|
|
Eigen::MatrixXd pHN;
|
|
Eigen::MatrixXd tHN;
|
|
const auto update = [&]()
|
|
{
|
|
// why does windows need this Eigen::VectorXd(...) ?
|
|
pHN = (pL*V).array().colwise() / Eigen::VectorXd(pM.diagonal()).array();
|
|
tHN = (tL*V).array().colwise() / Eigen::VectorXd(tM.diagonal()).array();
|
|
pHN *= 1.0/pHN.rowwise().norm().maxCoeff();
|
|
tHN *= 1.0/tHN.rowwise().norm().maxCoeff();
|
|
const auto was_face_based = vr.data_list[0].face_based;
|
|
Eigen::MatrixXd QV(V.rows()*2,3);
|
|
QV.topRows(V.rows()) = V;
|
|
if(use_poly)
|
|
{
|
|
printf("using polygon Laplacian\n");
|
|
QV.bottomRows(V.rows()) = V-pHN;
|
|
}else
|
|
{
|
|
printf("using triangle Laplacian\n");
|
|
QV.bottomRows(V.rows()) = V-tHN;
|
|
}
|
|
Eigen::MatrixXi QE(V.rows(),2);
|
|
for(int i = 0;i<V.rows();i++){ QE(i,0)=i;QE(i,1)=i+V.rows();}
|
|
vr.data_list[1].set_edges(QV,QE,Eigen::RowVector3d(1,1,1));
|
|
|
|
if(use_poly)
|
|
{
|
|
vr.data_list[0].show_lines = false;
|
|
if(show_edges)
|
|
{
|
|
vr.data_list[0].set_edges(V,E,Eigen::RowVector3d(0,0,0));
|
|
}else
|
|
{
|
|
vr.data_list[0].clear_edges();
|
|
}
|
|
}else
|
|
{
|
|
vr.data_list[0].clear_edges();
|
|
vr.data_list[0].show_lines = show_edges;
|
|
}
|
|
vr.data_list[0].set_face_based(was_face_based);
|
|
};
|
|
const double original_area = pM.diagonal().sum();
|
|
const auto recompute_M = [&]()
|
|
{
|
|
Eigen::SparseMatrix<double> _1,_2;
|
|
igl::cotmatrix(V,I,C,_1,pM,_2);
|
|
igl::massmatrix(V,F,igl::MASSMATRIX_TYPE_DEFAULT,tM);
|
|
V *= sqrt(original_area / pM.diagonal().sum());
|
|
};
|
|
const auto cmcf_step = [&]()
|
|
{
|
|
const Eigen::SparseMatrix<double> S =
|
|
use_poly? ((pM) - 0.05*(pL)): ((tM) - 0.05*(tL));
|
|
const Eigen::MatrixXd rhs = use_poly? pM*V : tM*V;
|
|
Eigen::SimplicialLLT<Eigen::SparseMatrix<double > > solver(S);
|
|
assert(solver.info() == Eigen::Success);
|
|
V = solver.solve(rhs).eval();
|
|
// recompute just mass matrices
|
|
recompute_M();
|
|
// center
|
|
V.rowwise() -= V.colwise().mean();
|
|
vr.data_list[0].set_vertices(V);
|
|
vr.data_list[0].compute_normals();
|
|
update();
|
|
};
|
|
vr.callback_key_pressed = [&](decltype(vr) &,unsigned int key, int mod)
|
|
{
|
|
switch(key)
|
|
{
|
|
case ' ': cmcf_step(); return true;
|
|
case 'R': case 'r': V=OV;recompute_M();vr.data_list[0].set_vertices(V);vr.data_list[0].compute_normals(); update();return true;
|
|
case 'P': case 'p': use_poly=!use_poly; update();return true;
|
|
case 'L': case 'l': show_edges=!show_edges; update();return true;
|
|
}
|
|
return false;
|
|
};
|
|
update();
|
|
vr.launch();
|
|
}
|
|
|