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.
105 lines
2.9 KiB
105 lines
2.9 KiB
1 year ago
|
#include <igl/barycenter.h>
|
||
|
|
#include <igl/cotmatrix.h>
|
||
|
|
#include <igl/doublearea.h>
|
||
|
|
#include <igl/grad.h>
|
||
|
|
#include <igl/jet.h>
|
||
|
|
#include <igl/massmatrix.h>
|
||
|
|
#include <igl/per_vertex_normals.h>
|
||
|
|
#include <igl/readDMAT.h>
|
||
|
|
#include <igl/readOFF.h>
|
||
|
|
#include <igl/repdiag.h>
|
||
|
|
#include <igl/opengl/glfw/Viewer.h>
|
||
|
|
|
||
|
|
#include <iostream>
|
||
|
|
|
||
|
|
Eigen::MatrixXd V,U;
|
||
|
|
Eigen::MatrixXi F;
|
||
|
|
Eigen::SparseMatrix<double> L;
|
||
|
|
igl::opengl::glfw::Viewer viewer;
|
||
|
|
|
||
|
|
int main(int argc, char *argv[])
|
||
|
|
{
|
||
|
|
using namespace Eigen;
|
||
|
|
using namespace std;
|
||
|
|
|
||
|
|
// Load a mesh in OFF format
|
||
|
|
igl::readOFF(TUTORIAL_SHARED_PATH "/cow.off", V, F);
|
||
|
|
|
||
|
|
// Compute Laplace-Beltrami operator: #V by #V
|
||
|
|
igl::cotmatrix(V,F,L);
|
||
|
|
|
||
|
|
// Alternative construction of same Laplacian
|
||
|
|
SparseMatrix<double> G,K;
|
||
|
|
// Gradient/Divergence
|
||
|
|
igl::grad(V,F,G);
|
||
|
|
// Diagonal per-triangle "mass matrix"
|
||
|
|
VectorXd dblA;
|
||
|
|
igl::doublearea(V,F,dblA);
|
||
|
|
// Place areas along diagonal #dim times
|
||
|
|
const auto & T = 1.*(dblA.replicate(3,1)*0.5).asDiagonal();
|
||
|
|
// Laplacian K built as discrete divergence of gradient or equivalently
|
||
|
|
// discrete Dirichelet energy Hessian
|
||
|
|
K = -G.transpose() * T * G;
|
||
|
|
cout<<"|K-L|: "<<(K-L).norm()<<endl;
|
||
|
|
|
||
|
|
const auto &key_down = [](igl::opengl::glfw::Viewer &viewer,unsigned char key,int mod)->bool
|
||
|
|
{
|
||
|
|
switch(key)
|
||
|
|
{
|
||
|
|
case 'r':
|
||
|
|
case 'R':
|
||
|
|
U = V;
|
||
|
|
break;
|
||
|
|
case ' ':
|
||
|
|
{
|
||
|
|
// Recompute just mass matrix on each step
|
||
|
|
SparseMatrix<double> M;
|
||
|
|
igl::massmatrix(U,F,igl::MASSMATRIX_TYPE_BARYCENTRIC,M);
|
||
|
|
// Solve (M-delta*L) U = M*U
|
||
|
|
const auto & S = (M - 0.001*L);
|
||
|
|
Eigen::SimplicialLLT<Eigen::SparseMatrix<double > > solver(S);
|
||
|
|
assert(solver.info() == Eigen::Success);
|
||
|
|
U = solver.solve(M*U).eval();
|
||
|
|
// Compute centroid and subtract (also important for numerics)
|
||
|
|
VectorXd dblA;
|
||
|
|
igl::doublearea(U,F,dblA);
|
||
|
|
double area = 0.5*dblA.sum();
|
||
|
|
MatrixXd BC;
|
||
|
|
igl::barycenter(U,F,BC);
|
||
|
|
RowVector3d centroid(0,0,0);
|
||
|
|
for(int i = 0;i<BC.rows();i++)
|
||
|
|
{
|
||
|
|
centroid += 0.5*dblA(i)/area*BC.row(i);
|
||
|
|
}
|
||
|
|
U.rowwise() -= centroid;
|
||
|
|
// Normalize to unit surface area (important for numerics)
|
||
|
|
U.array() /= sqrt(area);
|
||
|
|
break;
|
||
|
|
}
|
||
|
|
default:
|
||
|
|
return false;
|
||
|
|
}
|
||
|
|
// Send new positions, update normals, recenter
|
||
|
|
viewer.data().set_vertices(U);
|
||
|
|
viewer.data().compute_normals();
|
||
|
|
viewer.core().align_camera_center(U,F);
|
||
|
|
return true;
|
||
|
|
};
|
||
|
|
|
||
|
|
|
||
|
|
// Use original normals as pseudo-colors
|
||
|
|
MatrixXd N;
|
||
|
|
igl::per_vertex_normals(V,F,N);
|
||
|
|
MatrixXd C = N.rowwise().normalized().array()*0.5+0.5;
|
||
|
|
|
||
|
|
// Initialize smoothing with base mesh
|
||
|
|
U = V;
|
||
|
|
viewer.data().set_mesh(U, F);
|
||
|
|
viewer.data().set_colors(C);
|
||
|
|
viewer.callback_key_down = key_down;
|
||
|
|
|
||
|
|
cout<<"Press [space] to smooth."<<endl;;
|
||
|
|
cout<<"Press [r] to reset."<<endl;;
|
||
|
|
return viewer.launch();
|
||
|
|
}
|