First version of the implicit smoother
This commit is contained in:
parent
4dc3c08004
commit
d03a5a97d2
|
|
@ -0,0 +1,304 @@
|
|||
/****************************************************************************
|
||||
* VCGLib o o *
|
||||
* Visual and Computer Graphics Library o o *
|
||||
* _ O _ *
|
||||
* Copyright(C) 2004 \/)\/ *
|
||||
* Visual Computing Lab /\/| *
|
||||
* ISTI - Italian National Research Council | *
|
||||
* \ *
|
||||
* All rights reserved. *
|
||||
* *
|
||||
* This program is free software; you can redistribute it and/or modify *
|
||||
* it under the terms of the GNU General Public License as published by *
|
||||
* the Free Software Foundation; either version 2 of the License, or *
|
||||
* (at your option) any later version. *
|
||||
* *
|
||||
* This program is distributed in the hope that it will be useful, *
|
||||
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
|
||||
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
|
||||
* GNU General Public License (http://www.gnu.org/licenses/gpl.txt) *
|
||||
* for more details. *
|
||||
* *
|
||||
****************************************************************************/
|
||||
#ifndef __VCG_IMPLICIT_SMOOTHER
|
||||
#define __VCG_IMPLICIT_SMOOTHER
|
||||
|
||||
#include <eigenlib/Eigen/Sparse>
|
||||
#include <vcg/complex/algorithms/mesh_to_matrix.h>
|
||||
#include <vcg/complex/algorithms/update/quality.h>
|
||||
#include <vcg/complex/algorithms/smooth.h>
|
||||
|
||||
#define PENALTY 10000
|
||||
|
||||
namespace vcg{
|
||||
|
||||
|
||||
template <class MeshType>
|
||||
class ImplicitSmoother
|
||||
{
|
||||
typedef typename MeshType::FaceType FaceType;
|
||||
typedef typename MeshType::VertexType VertexType;
|
||||
typedef typename MeshType::CoordType CoordType;
|
||||
typedef typename MeshType::ScalarType ScalarType;
|
||||
typedef typename Eigen::Matrix<ScalarType, Eigen::Dynamic, Eigen::Dynamic> MatrixXm;
|
||||
|
||||
MeshType &to_smooth_mesh;
|
||||
|
||||
public:
|
||||
|
||||
struct FaceConstraint
|
||||
{
|
||||
int numF;
|
||||
std::vector<std::vector<ScalarType> > BarycentricW;
|
||||
CoordType TargetPos;
|
||||
ScalarType Weight;
|
||||
|
||||
FaceConstraint()
|
||||
{
|
||||
numF=-1;
|
||||
Weight=0;
|
||||
}
|
||||
|
||||
FaceConstraint(int _numF,const std::vector<std::vector<ScalarType> > &_BarycentricW,
|
||||
const CoordType &_TargetPos,ScalarType _Weight)
|
||||
{
|
||||
numF=_numF;
|
||||
BarycentricW= std::vector<std::vector<ScalarType> > (_BarycentricW.begin(),_BarycentricW.end());
|
||||
TargetPos=_TargetPos;
|
||||
Weight=_Weight;
|
||||
}
|
||||
};
|
||||
|
||||
struct SmoothParam
|
||||
{
|
||||
//the amount of smoothness, useful only if we set the mass matrix
|
||||
ScalarType lambda;
|
||||
//the use of mass matrix to keep the mesh close to its original position
|
||||
//(weighted per area distributed on vertices)
|
||||
bool useMassMatrix;
|
||||
//this bool is used to fix the border vertices of the mesh or not
|
||||
bool fixBorder;
|
||||
//this bool is used to set if cotangent weight is used, this flag to false means uniform laplacian
|
||||
bool useCotWeight;
|
||||
//the set of fixed vertices
|
||||
std::vector<int> FixedV;
|
||||
//the set of faces for barycentric constraints
|
||||
std::vector<FaceConstraint> ConstrainedF;
|
||||
|
||||
SmoothParam()
|
||||
{
|
||||
lambda=0.2;
|
||||
useMassMatrix=true;
|
||||
fixBorder=false;
|
||||
useCotWeight=false;
|
||||
}
|
||||
};
|
||||
|
||||
private:
|
||||
|
||||
// void GetMassMatrix(MeshType &mesh,
|
||||
// std::vector<std::pair<int,int> > &index,
|
||||
// std::vector<ScalarType> &entry)
|
||||
// {
|
||||
// entry.clear();
|
||||
// index.clear();
|
||||
|
||||
// //calculate area
|
||||
// vcg::tri::UpdateQuality<MeshType>::FaceArea(mesh);
|
||||
// //then distribute per vertex
|
||||
// vcg::tri::UpdateQuality<MeshType>::VertexFromFace(mesh);
|
||||
|
||||
// //3, one per coordinate
|
||||
// index.resize(mesh.vert.size()*3,-1);
|
||||
// entry.resize(mesh.vert.size()*3,0);
|
||||
|
||||
// //store the index and the scalar for the sparse matrix
|
||||
// for (size_t i=0;i<mesh.vert.size();i++)
|
||||
// {
|
||||
// for (size_t j=0;j<3;j++)
|
||||
// {
|
||||
// int currI=(i*3)+j;
|
||||
// index[currI]=currI;
|
||||
// entry[currI]=mesh.vert[i].Q();
|
||||
// }
|
||||
// }
|
||||
// }
|
||||
|
||||
void InitSparse(const std::vector<std::pair<int,int> > &Index,
|
||||
const std::vector<ScalarType> &Values,
|
||||
const size_t m,
|
||||
const size_t n,
|
||||
Eigen::SparseMatrix<ScalarType>& X)
|
||||
{
|
||||
|
||||
assert(Index.size()==Values.size());
|
||||
|
||||
std::vector<Eigen::Triplet<ScalarType> > IJV;
|
||||
IJV.reserve(Index.size());
|
||||
|
||||
for(size_t i= 0;i<Index.size();i++)
|
||||
{
|
||||
int row=Index[i].first;
|
||||
int col=Index[i].second;
|
||||
ScalarType val=Values[i];
|
||||
|
||||
assert(row<m);
|
||||
assert(col<n);
|
||||
|
||||
IJV.push_back(Eigen::Triplet<ScalarType>(row,col,val));
|
||||
}
|
||||
X.resize(m,n);
|
||||
X.setFromTriplets(IJV.begin(),IJV.end());
|
||||
}
|
||||
|
||||
int SystemSize(SmoothParam & SParam)
|
||||
{
|
||||
int basic_size=to_smooth_mesh.vert.size();
|
||||
int constr_size=SParam.ConstrainedF.size();
|
||||
return (basic_size+constr_size);
|
||||
}
|
||||
|
||||
void CollectHardConstraints(const SmoothParam &SParam,
|
||||
std::vector<std::pair<int,int> > &IndexC,
|
||||
std::vector<ScalarType> &WeightC)
|
||||
{
|
||||
std::vector<int> To_Fix;
|
||||
|
||||
//collect fixed vert
|
||||
if (SParam.fixBorder)
|
||||
{
|
||||
//add penalization constra
|
||||
for (int i=0;i<to_smooth_mesh.vert.size();i++)
|
||||
{
|
||||
if (!to_smooth_mesh.vert[i].IsB())continue;
|
||||
To_Fix.push_back(i);
|
||||
}
|
||||
}
|
||||
//add additional fixed vertices constraint
|
||||
To_Fix.insert(To_Fix.end(),SParam.FixedV.begin(),SParam.FixedV.end());
|
||||
|
||||
//sort and make them unique
|
||||
std::sort(To_Fix.begin(),To_Fix.end());
|
||||
typename std::vector<int>::iterator it= std::unique (To_Fix.begin(), To_Fix.end());
|
||||
To_Fix.resize( std::distance(To_Fix.begin(),it) );
|
||||
|
||||
for (size_t i=0;i<To_Fix.size();i++)
|
||||
{
|
||||
for (int j=0;j<3;j++)
|
||||
{
|
||||
int IndexV=(i*3)+j;
|
||||
IndexC.push_back(std::pair<int,int>(IndexV,IndexV));
|
||||
WeightC.push_back((ScalarType)PENALTY);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
public:
|
||||
|
||||
|
||||
|
||||
//static void Smooth
|
||||
// static void SmoothLIbIGL(MeshType &mesh,
|
||||
// ScalarType lambda=0.5)//0.0001)
|
||||
// {
|
||||
// Eigen::MatrixXi F;
|
||||
// Eigen::MatrixXd V,U;
|
||||
// vcg::tri::MeshToMatrix<MeshType>::GetTriMeshData(mesh,F,V);
|
||||
// U=V;
|
||||
|
||||
// Eigen::SparseMatrix<ScalarType> L,M;
|
||||
// igl::cotmatrix(V,F,L);
|
||||
|
||||
// //compute the mass matrix
|
||||
// static bool computed=false;
|
||||
// igl::massmatrix(V,F,igl::MASSMATRIX_TYPE_BARYCENTRIC,M);
|
||||
// M.setIdentity();
|
||||
// computed=true;
|
||||
|
||||
// //const auto & S = (M - 0.001*L);
|
||||
// Eigen::SparseMatrix<double> S = (M - lambda*L);
|
||||
|
||||
|
||||
// Eigen::SimplicialLLT<Eigen::SparseMatrix<double > > solver(S);
|
||||
// assert(solver.info() == Eigen::Success);
|
||||
|
||||
// U = solver.solve(M*U).eval();
|
||||
// // Normalize to unit surface area (important for numerics)
|
||||
// U.array() /= sqrt(M.diagonal().array().sum());
|
||||
|
||||
// //then assing per vertex positions
|
||||
// printf("smoothed \n");
|
||||
// fflush(stdout);
|
||||
// for (size_t i=0;i<mesh.vert.size();i++)
|
||||
// {
|
||||
// mesh.vert[i].P().X()=U(i,0);
|
||||
// mesh.vert[i].P().Y()=U(i,1);
|
||||
// mesh.vert[i].P().Z()=U(i,2);
|
||||
// }
|
||||
// }
|
||||
|
||||
void Smooth(SmoothParam &SParam)
|
||||
{
|
||||
|
||||
//the laplacian and the mass matrix
|
||||
Eigen::SparseMatrix<ScalarType> L,M;
|
||||
|
||||
//initialize the mass matrix
|
||||
std::vector<std::pair<int,int> > IndexM;
|
||||
std::vector<ScalarType> ValuesM;
|
||||
|
||||
//add the entries for mass matrix
|
||||
if (SParam.useMassMatrix) MeshToMatrix<MeshType>::MassMatrixEntry(to_smooth_mesh,IndexM,ValuesM);
|
||||
|
||||
//then also collect hard constraints
|
||||
//CollectHardConstraints(SParam,IndexM,ValuesM);
|
||||
|
||||
//initialize sparse mass matrix
|
||||
InitSparse(IndexM,ValuesM,to_smooth_mesh.vert.size()*3,to_smooth_mesh.vert.size()*3,M);
|
||||
|
||||
//get the entries for laplacian matrix
|
||||
std::vector<std::pair<int,int> > IndexL;
|
||||
std::vector<ScalarType> ValuesL;
|
||||
MeshToMatrix<MeshType>::GetLaplacianMatrix(to_smooth_mesh,IndexL,ValuesL,false);//SParam.useCotWeight);
|
||||
|
||||
//initialize sparse laplacian matrix
|
||||
InitSparse(IndexL,ValuesL,to_smooth_mesh.vert.size()*3,to_smooth_mesh.vert.size()*3,L);
|
||||
|
||||
//then solve the system
|
||||
Eigen::SparseMatrix<ScalarType> S = (M + SParam.lambda*L);
|
||||
|
||||
//SimplicialLDLT
|
||||
Eigen::SimplicialLDLT<Eigen::SparseMatrix<ScalarType > > solver(S);
|
||||
printf("output %d \n",solver.info());
|
||||
fflush(stdout);
|
||||
assert(solver.info() == Eigen::Success);
|
||||
|
||||
int size=to_smooth_mesh.vert.size()*3;
|
||||
MatrixXm V(size,1);
|
||||
|
||||
for (size_t i=0;i<to_smooth_mesh.vert.size();i++)
|
||||
{
|
||||
int index=i*3;
|
||||
assert(index<size);
|
||||
V(index,0)=to_smooth_mesh.vert[i].P().X();
|
||||
V(index+1,0)=to_smooth_mesh.vert[i].P().Y();
|
||||
V(index+2,0)=to_smooth_mesh.vert[i].P().Z();
|
||||
}
|
||||
|
||||
V = solver.solve(M*V).eval();
|
||||
|
||||
for (size_t i=0;i<to_smooth_mesh.vert.size();i++)
|
||||
{
|
||||
int index=i*3;
|
||||
to_smooth_mesh.vert[i].P().X()=V(index,0);
|
||||
to_smooth_mesh.vert[i].P().Y()=V(index+1,0);
|
||||
to_smooth_mesh.vert[i].P().Z()=V(index+2,0);
|
||||
}
|
||||
}
|
||||
|
||||
ImplicitSmoother(MeshType &_to_smooth_mesh):to_smooth_mesh(_to_smooth_mesh){}
|
||||
};
|
||||
|
||||
}//end namespace vcg
|
||||
|
||||
#endif
|
||||
|
|
@ -25,23 +25,26 @@
|
|||
|
||||
#include <vcg/complex/complex.h>
|
||||
#include <vcg/complex/algorithms/update/topology.h>
|
||||
#include <vcg/complex/algorithms/update/quality.h>
|
||||
#include <vcg/complex/algorithms/harmonic.h>
|
||||
|
||||
using namespace std;
|
||||
|
||||
namespace vcg {
|
||||
namespace tri {
|
||||
template < typename TriMeshType >
|
||||
template < typename MeshType >
|
||||
class MeshToMatrix
|
||||
{
|
||||
|
||||
// define types
|
||||
|
||||
typedef typename TriMeshType::FaceType FaceType;
|
||||
typedef typename TriMeshType::VertexType VertexType;
|
||||
typedef typename TriMeshType::CoordType CoordType;
|
||||
typedef typename TriMeshType::ScalarType ScalarType;
|
||||
typedef typename MeshType::FaceType FaceType;
|
||||
typedef typename MeshType::VertexType VertexType;
|
||||
typedef typename MeshType::CoordType CoordType;
|
||||
typedef typename MeshType::ScalarType ScalarType;
|
||||
typedef typename Eigen::Matrix<ScalarType, Eigen::Dynamic, Eigen::Dynamic> MatrixXm;
|
||||
|
||||
static void GetTriEdgeAdjacency(const Eigen::MatrixXd& V,
|
||||
static void GetTriEdgeAdjacency(const MatrixXm& V,
|
||||
const Eigen::MatrixXi& F,
|
||||
Eigen::MatrixXi& EV,
|
||||
Eigen::MatrixXi& FE,
|
||||
|
|
@ -128,13 +131,13 @@ class MeshToMatrix
|
|||
public:
|
||||
|
||||
// return mesh as vector of vertices and faces
|
||||
static void GetTriMeshData(const TriMeshType &mesh,
|
||||
static void GetTriMeshData(const MeshType &mesh,
|
||||
Eigen::MatrixXi &faces,
|
||||
Eigen::MatrixXd &vert)
|
||||
MatrixXm &vert)
|
||||
{
|
||||
tri::RequireCompactness(mesh);
|
||||
// create eigen matrix of vertices
|
||||
vert=Eigen::MatrixXd(mesh.VN(), 3);
|
||||
vert=MatrixXm(mesh.VN(), 3);
|
||||
|
||||
// copy vertices
|
||||
for (int i = 0; i < mesh.VN(); i++)
|
||||
|
|
@ -151,13 +154,13 @@ public:
|
|||
}
|
||||
|
||||
// return normals of the mesh
|
||||
static void GetNormalData(const TriMeshType &mesh,
|
||||
Eigen::MatrixXd &Nvert,
|
||||
Eigen::MatrixXd &Nface)
|
||||
static void GetNormalData(const MeshType &mesh,
|
||||
MatrixXm &Nvert,
|
||||
MatrixXm &Nface)
|
||||
{
|
||||
// create eigen matrix of vertices
|
||||
Nvert=Eigen::MatrixXd(mesh.VN(), 3);
|
||||
Nface=Eigen::MatrixXd(mesh.FN(), 3);
|
||||
Nvert=MatrixXm(mesh.VN(), 3);
|
||||
Nface=MatrixXm(mesh.FN(), 3);
|
||||
|
||||
// per vertices normals
|
||||
for (int i = 0; i < mesh.VN(); i++)
|
||||
|
|
@ -171,11 +174,11 @@ public:
|
|||
}
|
||||
|
||||
// get face to face adjacency
|
||||
static void GetTriFFAdjacency(TriMeshType &mesh,
|
||||
static void GetTriFFAdjacency(MeshType &mesh,
|
||||
Eigen::MatrixXi &FFp,
|
||||
Eigen::MatrixXi &FFi)
|
||||
{
|
||||
tri::UpdateTopology<TriMeshType>::FaceFace(mesh);
|
||||
tri::UpdateTopology<MeshType>::FaceFace(mesh);
|
||||
FFp = Eigen::MatrixXi(mesh.FN(),3);
|
||||
FFi = Eigen::MatrixXi(mesh.FN(),3);
|
||||
|
||||
|
|
@ -197,13 +200,13 @@ public:
|
|||
}
|
||||
|
||||
// get edge to face and edge to vertex adjacency
|
||||
static void GetTriEdgeAdjacency(const TriMeshType &mesh,
|
||||
static void GetTriEdgeAdjacency(const MeshType &mesh,
|
||||
Eigen::MatrixXi& EV,
|
||||
Eigen::MatrixXi& FE,
|
||||
Eigen::MatrixXi& EF)
|
||||
{
|
||||
Eigen::MatrixXi faces;
|
||||
Eigen::MatrixXd vert;
|
||||
MatrixXm vert;
|
||||
GetTriMeshData(mesh,faces,vert);
|
||||
GetTriEdgeAdjacency(vert,faces,EV,FE,EF);
|
||||
}
|
||||
|
|
@ -213,6 +216,85 @@ public:
|
|||
Eigen::Vector3d ret(v[0],v[1],v[2]);
|
||||
return ret;
|
||||
}
|
||||
|
||||
|
||||
static void MassMatrixEntry(MeshType &mesh,
|
||||
std::vector<std::pair<int,int> > &index,
|
||||
std::vector<ScalarType> &entry)
|
||||
{
|
||||
//calculate area
|
||||
UpdateQuality<MeshType>::FaceArea(mesh);
|
||||
//then distribute per vertex
|
||||
UpdateQuality<MeshType>::VertexFromFace(mesh);
|
||||
|
||||
std::pair<float,float> minmax=Stat<MeshType>::ComputePerVertexQualityMinMax(mesh);
|
||||
|
||||
//store the index and the scalar for the sparse matrix
|
||||
for (size_t i=0;i<mesh.vert.size();i++)
|
||||
{
|
||||
for (size_t j=0;j<3;j++)
|
||||
{
|
||||
int currI=(i*3)+j;
|
||||
index.push_back(std::pair<int,int>(currI,currI));
|
||||
entry.push_back(mesh.vert[i].Q()/(ScalarType)minmax.second);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
static void GetLaplacianEntry(MeshType &mesh,
|
||||
FaceType &f,
|
||||
std::vector<std::pair<int,int> > &index,
|
||||
std::vector<ScalarType> &entry,
|
||||
bool cotangent)
|
||||
{
|
||||
for (int i=0;i<3;i++)
|
||||
{
|
||||
|
||||
ScalarType weight = 1;
|
||||
if (cotangent)
|
||||
{
|
||||
weight=Harmonic<MeshType>::template CotangentWeight<ScalarType>(f,i);
|
||||
}
|
||||
|
||||
//get the index of the vertices
|
||||
int indexV0=Index(mesh,f.V0(i));
|
||||
int indexV1=Index(mesh,f.V1(i));
|
||||
|
||||
//then assemble the matrix
|
||||
for (int j=0;j<3;j++)
|
||||
{
|
||||
//multiply by 3 and add the component
|
||||
int currI0=(indexV0*3)+j;
|
||||
int currI1=(indexV1*3)+j;
|
||||
|
||||
index.push_back(std::pair<int,int>(currI0,currI0));
|
||||
entry.push_back(weight);
|
||||
index.push_back(std::pair<int,int>(currI0,currI1));
|
||||
entry.push_back(-weight);
|
||||
|
||||
index.push_back(std::pair<int,int>(currI1,currI1));
|
||||
entry.push_back(weight);
|
||||
index.push_back(std::pair<int,int>(currI1,currI0));
|
||||
entry.push_back(-weight);
|
||||
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
static void GetLaplacianMatrix(MeshType &mesh,
|
||||
std::vector<std::pair<int,int> > &index,
|
||||
std::vector<ScalarType> &entry,
|
||||
bool cotangent)
|
||||
{
|
||||
//store the index and the scalar for the sparse matrix
|
||||
for (size_t i=0;i<mesh.face.size();i++)
|
||||
GetLaplacianEntry(mesh,mesh.face[i],index,entry,cotangent);
|
||||
}
|
||||
|
||||
|
||||
|
||||
};
|
||||
|
||||
} // end namespace tri
|
||||
|
|
|
|||
Loading…
Reference in New Issue