mirror of
https://github.com/lucaspalomodevelop/meshlab.git
synced 2026-03-13 00:07:24 +00:00
version 1.1
Added matrix caching, major speedup in cotan weights computation
This commit is contained in:
davidboening
parent
a3879881cc
commit
f8f213ea4e
@ -173,86 +173,8 @@ std::map<std::string, QVariant> FilterHeatGeodesicPlugin::applyFilter(const QAct
|
||||
mm.updateDataMask(MeshModel::MM_VERTQUALITY);
|
||||
mm.updateDataMask(MeshModel::MM_VERTCOLOR);
|
||||
CMeshO &mesh = mm.cm;
|
||||
|
||||
cb(100*0/10, "Updating Topology and Computing Face Normals...");
|
||||
vcg::tri::UpdateTopology<CMeshO>::VertexFace(mesh);
|
||||
vcg::tri::UpdateTopology<CMeshO>::FaceFace(mesh);
|
||||
vcg::tri::UpdateNormal<CMeshO>::PerFaceNormalized(mesh);
|
||||
|
||||
cb(100*1/10, "Computing Boundary Conditions...");
|
||||
Eigen::VectorXd initialConditions(mesh.VN());
|
||||
int selection_count = 0;
|
||||
for(int i=0; i < mesh.vert.size(); i++){
|
||||
if (mesh.vert[i].IsS()){
|
||||
initialConditions[i] = 1;
|
||||
++selection_count;
|
||||
}
|
||||
else {
|
||||
initialConditions[i] = 0;
|
||||
}
|
||||
}
|
||||
if (selection_count < 1){
|
||||
log("Warning: no vertices are selected! aborting computation.");
|
||||
break;
|
||||
}
|
||||
|
||||
|
||||
cb(100*2/10, "Building Linear System (Heatflow)...");
|
||||
Eigen::SparseMatrix<double> mass(mesh.VN(), mesh.VN());
|
||||
buildMassMatrix(mesh, mass);
|
||||
|
||||
Eigen::SparseMatrix<double> cotanOperator(mesh.VN(), mesh.VN());
|
||||
buildCotanMatrix(mesh, cotanOperator);
|
||||
|
||||
double avg_edge_len = computeAverageEdgeLength(mesh);
|
||||
double timestep = parameters.getFloat("m") * avg_edge_len * avg_edge_len;
|
||||
Eigen::SparseMatrix<double> system1(mesh.VN(), mesh.VN());
|
||||
system1 = mass - timestep * cotanOperator;
|
||||
|
||||
Eigen::SimplicialLDLT<Eigen::SparseMatrix<double>> solver;
|
||||
|
||||
cb(100*3/10, "Computing Matrix Factorization (Heatflow)...");
|
||||
solver.compute(system1);
|
||||
if(solver.info() != Eigen::Success) {
|
||||
log("Error: Factorization Failed (Heatflow).");
|
||||
}
|
||||
cb(100*4/10, "Solving Linear System (Heatflow)...");
|
||||
Eigen::VectorXd heatflow = solver.solve(initialConditions); // (VN)
|
||||
if(solver.info() != Eigen::Success) {
|
||||
log("Error: Solver Failed (Heatflow).");
|
||||
}
|
||||
cb(100*5/10, "Computing Heat Gradient...");
|
||||
Eigen::MatrixX3d heatGradient = computeVertexGradient(mesh, heatflow); // (FN, 3)
|
||||
|
||||
Eigen::MatrixX3d normalizedVectorField = normalizeVectorField(-heatGradient); // (FN, 3)
|
||||
cb(100*6/10, "Computing Divergence...");
|
||||
Eigen::VectorXd divergence = computeVertexDivergence(mesh, normalizedVectorField); // (VN)
|
||||
|
||||
cb(100*7/10, "Building Linear System (Geodesic)...");
|
||||
Eigen::SparseMatrix<double> system2(mesh.VN(), mesh.VN());
|
||||
system2 = cotanOperator; //+ 1e-6 * Eigen::Matrix<double,-1,-1>::Identity(mesh.VN(), mesh.VN());
|
||||
|
||||
cb(100*8/10, "Computing Matrix Factorization (Geodesic)...");
|
||||
solver.compute(system2);
|
||||
if(solver.info() != Eigen::Success) {
|
||||
log("Error: Factorization Failed (Geodesic).");
|
||||
}
|
||||
cb(100*9/10, "Solving Linear System (Geodesic)...");
|
||||
Eigen::VectorXd geodesicDistance = solver.solve(divergence);
|
||||
if(solver.info() != Eigen::Success) {
|
||||
log("Error: Solver Failed (Geodesic).");
|
||||
}
|
||||
cb(100*10/10, "Saving Geodesic Distance...");
|
||||
// invert and shift
|
||||
geodesicDistance.array() *= -1; // no clue as to why this needs to be here
|
||||
geodesicDistance.array() -= geodesicDistance.minCoeff();
|
||||
|
||||
// set geodesic distance as quality
|
||||
for(int i=0; i < mesh.vert.size(); i++){
|
||||
mesh.vert[i].Q() = geodesicDistance(i);
|
||||
}
|
||||
|
||||
vcg::tri::UpdateColor<CMeshO>::PerVertexQualityRamp(mesh);
|
||||
// TODO: compact all vectors
|
||||
computeHeatGeodesicFromSelection(mesh, cb, parameters.getFloat("m"));
|
||||
}
|
||||
break;
|
||||
default :
|
||||
@ -261,4 +183,110 @@ std::map<std::string, QVariant> FilterHeatGeodesicPlugin::applyFilter(const QAct
|
||||
return std::map<std::string, QVariant>();
|
||||
}
|
||||
|
||||
inline void FilterHeatGeodesicPlugin::computeHeatGeodesicFromSelection(CMeshO& mesh, vcg::CallBackPos* cb, float m){
|
||||
// build boundary conditions
|
||||
Eigen::VectorXd boundaryConditions(mesh.VN());
|
||||
int selection_count = 0;
|
||||
for(int i=0; i < mesh.VN(); i++){
|
||||
boundaryConditions(i) = mesh.vert[i].IsS() ? (++selection_count, 1) : 0;
|
||||
}
|
||||
if (selection_count < 1){
|
||||
log("Warning: no vertices are selected! aborting computation.");
|
||||
return;
|
||||
}
|
||||
|
||||
// update topology and face normals
|
||||
vcg::tri::UpdateTopology<CMeshO>::VertexFace(mesh);
|
||||
vcg::tri::UpdateTopology<CMeshO>::FaceFace(mesh);
|
||||
vcg::tri::UpdateNormal<CMeshO>::PerFaceNormalized(mesh);
|
||||
|
||||
// state variables
|
||||
Eigen::SparseMatrix<double> massMatrix;
|
||||
Eigen::SparseMatrix<double> cotanMatrix;
|
||||
double avg_edge_len;
|
||||
|
||||
typedef std::tuple<
|
||||
Eigen::SparseMatrix<double>, // system1 fact
|
||||
Eigen::SparseMatrix<double>, // system2 fact
|
||||
double // avg edge len
|
||||
> HeatMethodData;
|
||||
|
||||
// recover state if it exists
|
||||
if (!vcg::tri::HasPerMeshAttribute(mesh, "HeatMethodData")){
|
||||
// current solution saves matrices not factorizations
|
||||
// if we save factorization it is best to at least also save cotanMatrix
|
||||
// to avoid an expensive reconstruction when parameter("m") changes
|
||||
buildMassMatrix(mesh, massMatrix);
|
||||
// this step is very slow (95% of total time) hence we pass the callback
|
||||
// to update progress and avoid freezes
|
||||
buildCotanLowerTriMatrix(mesh, cotanMatrix, cb);
|
||||
avg_edge_len = computeAverageEdgeLength(mesh);
|
||||
|
||||
HeatMethodData &HMData = vcg::tri::Allocator<CMeshO>::GetPerMeshAttribute<HeatMethodData>(mesh, std::string("HeatMethodData"))();
|
||||
HMData = HeatMethodData(massMatrix, cotanMatrix, avg_edge_len);
|
||||
}
|
||||
else {
|
||||
HeatMethodData &HMData = vcg::tri::Allocator<CMeshO>::GetPerMeshAttribute<HeatMethodData>(mesh, std::string("HeatMethodData"))();
|
||||
std::tie (massMatrix, cotanMatrix, avg_edge_len) = HMData;
|
||||
}
|
||||
|
||||
// cholesky type solver
|
||||
Eigen::SimplicialLDLT<Eigen::SparseMatrix<double>> solver;
|
||||
|
||||
// build system 1
|
||||
double timestep = m * avg_edge_len * avg_edge_len;
|
||||
Eigen::SparseMatrix<double> system1(mesh.VN(), mesh.VN());
|
||||
system1 = massMatrix - timestep * cotanMatrix;
|
||||
|
||||
cb(91, "Computing Factorization 1...");
|
||||
|
||||
// factorize system 1
|
||||
solver.compute(system1);
|
||||
if(solver.info() != Eigen::Success) {
|
||||
log("Error: Factorization Failed (Heatflow).");
|
||||
}
|
||||
|
||||
cb(93, "Solving System 1...");
|
||||
// solve system 1
|
||||
Eigen::VectorXd heatflow = solver.solve(boundaryConditions); // (VN)
|
||||
if(solver.info() != Eigen::Success) {
|
||||
log("Error: Solver Failed (Heatflow).");
|
||||
}
|
||||
|
||||
cb(95, "Computing Gradient, VectorField, Divergence...");
|
||||
// intermediate steps
|
||||
Eigen::MatrixX3d heatGradient = computeVertexGradient(mesh, heatflow); // (FN, 3)
|
||||
Eigen::MatrixX3d normalizedVectorField = normalizeVectorField(-heatGradient); // (FN, 3)
|
||||
Eigen::VectorXd divergence = computeVertexDivergence(mesh, normalizedVectorField); // (VN)
|
||||
|
||||
// build system 2
|
||||
Eigen::SparseMatrix<double> system2(mesh.VN(), mesh.VN());
|
||||
system2 = cotanMatrix; //+ 1e-6 * Eigen::Matrix<double,-1,-1>::Identity(mesh.VN(), mesh.VN());
|
||||
|
||||
cb(96, "Computing Factorization 2...");
|
||||
// factorize system 2
|
||||
solver.compute(system2);
|
||||
if(solver.info() != Eigen::Success) {
|
||||
log("Error: Factorization Failed (Geodesic).");
|
||||
}
|
||||
|
||||
cb(97, "Solving System 2...");
|
||||
// solve system 2
|
||||
Eigen::VectorXd geodesicDistance = solver.solve(divergence);
|
||||
|
||||
if(solver.info() != Eigen::Success) {
|
||||
log("Error: Solver Failed (Geodesic).");
|
||||
}
|
||||
|
||||
// shift to impose boundary conditions (dist(d) = 0 \forall d \in init_cond)
|
||||
geodesicDistance.array() -= geodesicDistance.minCoeff();
|
||||
|
||||
cb(99, "Updating Quality and Color...");
|
||||
// set geodesic distance as quality and color
|
||||
for(int i=0; i < mesh.vert.size(); i++){
|
||||
mesh.vert[i].Q() = geodesicDistance(i);
|
||||
}
|
||||
vcg::tri::UpdateColor<CMeshO>::PerVertexQualityRamp(mesh);
|
||||
}
|
||||
|
||||
MESHLAB_PLUGIN_NAME_EXPORTER(FilterSamplePlugin)
|
||||
|
||||
@ -33,7 +33,7 @@ public:
|
||||
vcg::CallBackPos * cb);
|
||||
|
||||
private:
|
||||
void computeHeatGeodesicFromSelection(MeshDocument& mesh, vcg::CallBackPos* cb, float m);
|
||||
void computeHeatGeodesicFromSelection(CMeshO& mesh, vcg::CallBackPos* cb, float m);
|
||||
};
|
||||
|
||||
#endif
|
||||
|
||||
@ -9,30 +9,26 @@
|
||||
#include <common/plugins/interfaces/filter_plugin.h>
|
||||
|
||||
|
||||
// foward declarations
|
||||
inline Eigen::Vector3d toEigen(const vcg::Point3f& p);
|
||||
inline double cotan(const Eigen::Vector3d& v0, const Eigen::Vector3d& v1);
|
||||
inline void buildMassMatrix(CMeshO &mesh, Eigen::SparseMatrix<double> &mass);
|
||||
inline void buildCotanMatrix(CMeshO &mesh, Eigen::SparseMatrix<double> &cotanOperator);
|
||||
inline double computeAverageEdgeLength(CMeshO &mesh);
|
||||
inline Eigen::MatrixX3d computeVertexGradient(CMeshO &mesh, const Eigen::VectorXd &heat);
|
||||
inline Eigen::VectorXd computeVertexDivergence(CMeshO &mesh, const Eigen::MatrixX3d &field);
|
||||
inline Eigen::MatrixX3d normalizeVectorField(const Eigen::MatrixX3d &field);
|
||||
// LOCAL FUNCTIONS
|
||||
|
||||
inline Eigen::Vector3d toEigen(const vcg::Point3f& p)
|
||||
{
|
||||
return Eigen::Vector3d(p.X(), p.Y(), p.Z());
|
||||
};
|
||||
namespace {
|
||||
inline Eigen::Vector3d toEigen(const vcg::Point3f& p)
|
||||
{
|
||||
return Eigen::Vector3d(p.X(), p.Y(), p.Z());
|
||||
};
|
||||
|
||||
|
||||
inline double cotan(const Eigen::Vector3d& v0, const Eigen::Vector3d& v1)
|
||||
{
|
||||
// cos(theta) / sin(theta)
|
||||
return v0.dot(v1) / v0.cross(v1).norm();
|
||||
};
|
||||
inline double cotan(const Eigen::Vector3d& v0, const Eigen::Vector3d& v1)
|
||||
{
|
||||
// cos(theta) / sin(theta)
|
||||
return v0.dot(v1) / v0.cross(v1).norm();
|
||||
};
|
||||
}
|
||||
|
||||
// GLOBAL FUNCTIONS
|
||||
|
||||
inline void buildMassMatrix(CMeshO &mesh, Eigen::SparseMatrix<double> &mass){
|
||||
mass.resize(mesh.VN(), mesh.VN());
|
||||
// compute area of all faces
|
||||
for (CMeshO::FaceIterator fi = mesh.face.begin(); fi != mesh.face.end(); ++fi)
|
||||
{
|
||||
@ -55,7 +51,7 @@ inline void buildMassMatrix(CMeshO &mesh, Eigen::SparseMatrix<double> &mass){
|
||||
std::vector<CMeshO::FaceType*> faces;
|
||||
std::vector<int> indices;
|
||||
vcg::face::VFStarVF<CMeshO::FaceType>(vp, faces, indices);
|
||||
|
||||
|
||||
double area = 0;
|
||||
for (int j = 0; j < faces.size(); ++j)
|
||||
{
|
||||
@ -67,55 +63,65 @@ inline void buildMassMatrix(CMeshO &mesh, Eigen::SparseMatrix<double> &mass){
|
||||
}
|
||||
|
||||
|
||||
inline void buildCotanMatrix(CMeshO &mesh, Eigen::SparseMatrix<double> &cotanOperator){
|
||||
inline void buildCotanLowerTriMatrix(CMeshO &mesh, Eigen::SparseMatrix<double> &cotanOperator, vcg::CallBackPos* cb = NULL){
|
||||
cotanOperator.resize(mesh.VN(), mesh.VN());
|
||||
// initialize a hashtable from vertex pointers to ids
|
||||
std::unordered_map<CMeshO::VertexType*, int> vertex_ids;
|
||||
for (int i = 0; i < mesh.VN(); ++i){
|
||||
vertex_ids[&mesh.vert[i]] = i;
|
||||
}
|
||||
|
||||
// iterate over all vertices to fill cotan matrix
|
||||
for (int i = 0; i < mesh.VN(); ++i){
|
||||
CMeshO::VertexType *vp = &mesh.vert[i];
|
||||
CMeshO::FaceType *fp = vp->VFp();
|
||||
vcg::face::Pos<CMeshO::FaceType> pos(fp, vp);
|
||||
vcg::face::Pos<CMeshO::FaceType> start(fp, vp);
|
||||
// iterate over all incident edges of vp
|
||||
do {
|
||||
// get the vertex opposite to vp
|
||||
pos.FlipV();
|
||||
CMeshO::VertexType *vo = pos.V();
|
||||
// move to left vertex
|
||||
pos.FlipE();pos.FlipV();
|
||||
CMeshO::VertexType *vl = pos.V();
|
||||
// move back then to right vertex
|
||||
pos.FlipV();pos.FlipE(); // back to vo
|
||||
pos.FlipF();pos.FlipE();pos.FlipV();
|
||||
CMeshO::VertexType *vr = pos.V();
|
||||
pos.FlipV();pos.FlipE();pos.FlipF();pos.FlipV(); // back to vp
|
||||
// progress bar variables
|
||||
int update_size = mesh.FN() / 90; // 90 updates (90% weight of the progress bar)
|
||||
int progress = 0;
|
||||
|
||||
// compute cotan of left edges and right edges
|
||||
Eigen::Vector3d elf = toEigen(vo->P() - vl->P()); // far left edge
|
||||
Eigen::Vector3d eln = toEigen(vp->P() - vl->P()); // near left edge
|
||||
Eigen::Vector3d erf = toEigen(vp->P() - vr->P()); // far right edge
|
||||
Eigen::Vector3d ern = toEigen(vo->P() - vr->P()); // near right edge
|
||||
// compute cotan weights
|
||||
for (int i = 0; i < mesh.FN(); ++i){
|
||||
// progress bar update
|
||||
if (cb != NULL && i % update_size == 89){
|
||||
cb(++progress, "Computing Cotan Weights...");
|
||||
}
|
||||
CMeshO::FaceType* fi = &mesh.face[i];
|
||||
|
||||
double cotan_l = cotan(elf, eln);
|
||||
double cotan_r = cotan(ern, erf);
|
||||
CMeshO::VertexType* v0 = fi->V(0);
|
||||
CMeshO::VertexType* v1 = fi->V(1);
|
||||
CMeshO::VertexType* v2 = fi->V(2);
|
||||
|
||||
// add to the matrix
|
||||
cotanOperator.coeffRef(vertex_ids[vp], vertex_ids[vo]) = (cotan_l + cotan_r)/2;
|
||||
vcg::Point3f p0 = v0->P();
|
||||
vcg::Point3f p1 = v1->P();
|
||||
vcg::Point3f p2 = v2->P();
|
||||
|
||||
// move to the next edge
|
||||
pos.FlipF();pos.FlipE();
|
||||
} while (pos != start);
|
||||
Eigen::Vector3d e0 = toEigen(p2 - p1);
|
||||
Eigen::Vector3d e1 = toEigen(p0 - p2);
|
||||
Eigen::Vector3d e2 = toEigen(p1 - p0);
|
||||
|
||||
// first edge is inverted to get correct orientation
|
||||
double alpha0 = cotan(-e1, e2) / 2;
|
||||
double alpha1 = cotan(-e2, e0) / 2;
|
||||
double alpha2 = cotan(-e0, e1) / 2;
|
||||
|
||||
int i0 = vertex_ids[v0];
|
||||
int i1 = vertex_ids[v1];
|
||||
int i2 = vertex_ids[v2];
|
||||
|
||||
// save only lower triangular part
|
||||
if (i0 > i1)
|
||||
cotanOperator.coeffRef(i0, i1) += alpha2;
|
||||
else
|
||||
cotanOperator.coeffRef(i1, i0) += alpha2;
|
||||
if (i0 > i2)
|
||||
cotanOperator.coeffRef(i0, i2) += alpha1;
|
||||
else
|
||||
cotanOperator.coeffRef(i2, i0) += alpha1;
|
||||
if (i1 > i2)
|
||||
cotanOperator.coeffRef(i1, i2) += alpha0;
|
||||
else
|
||||
cotanOperator.coeffRef(i2, i1) += alpha0;
|
||||
|
||||
cotanOperator.coeffRef(i0, i0) -= (alpha1 + alpha2);
|
||||
cotanOperator.coeffRef(i1, i1) -= (alpha0 + alpha2);
|
||||
cotanOperator.coeffRef(i2, i2) -= (alpha0 + alpha1);
|
||||
}
|
||||
|
||||
// compute diagonal entries
|
||||
for (int i = 0; i < mesh.VN(); ++i){
|
||||
cotanOperator.coeffRef(i, i) = -cotanOperator.row(i).sum();
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
|
||||
@ -156,8 +162,7 @@ inline Eigen::MatrixX3d computeVertexGradient(CMeshO &mesh, const Eigen::VectorX
|
||||
n /= n.norm();
|
||||
// face area
|
||||
double faceArea = fp->Q();
|
||||
// (ORDERING): edge unit vectors (assuming counter-clockwise ordering)
|
||||
// note if the ordering is clockwise, the gradient will point in the opposite direction
|
||||
// edge unit vectors (counter-clockwise)
|
||||
Eigen::Vector3d e0 = toEigen(p2 - p1);
|
||||
e0 /= e0.norm();
|
||||
Eigen::Vector3d e1 = toEigen(p0 - p2);
|
||||
@ -170,13 +175,13 @@ inline Eigen::MatrixX3d computeVertexGradient(CMeshO &mesh, const Eigen::VectorX
|
||||
Eigen::Vector3d g2 = n.cross(e2); //v2 grad
|
||||
|
||||
// add vertex gradient contributions
|
||||
Eigen::Vector3d total_grad = (
|
||||
g0 * heat(vertex_ids[fp->V(0)]) +
|
||||
g1 * heat(vertex_ids[fp->V(1)]) +
|
||||
g2 * heat(vertex_ids[fp->V(2)])
|
||||
Eigen::Vector3d tri_grad = (
|
||||
g0 * heat(vertex_ids[fp->V(0)]) +
|
||||
g1 * heat(vertex_ids[fp->V(1)]) +
|
||||
g2 * heat(vertex_ids[fp->V(2)])
|
||||
) / (2 * faceArea);
|
||||
|
||||
heatGradientField.row(i) = total_grad;
|
||||
heatGradientField.row(i) = tri_grad;
|
||||
}
|
||||
return heatGradientField;
|
||||
}
|
||||
@ -227,14 +232,14 @@ inline Eigen::VectorXd computeVertexDivergence(CMeshO &mesh, const Eigen::Matrix
|
||||
el = toEigen(p0 - p1); //-e2
|
||||
er = toEigen(p2 - p1); //e0
|
||||
eo = toEigen(p0 - p2); //e1
|
||||
} else if (index == 2){
|
||||
} else if (index == 2){
|
||||
el = toEigen(p1 - p2); //-e0
|
||||
er = toEigen(p0 - p2); //e1
|
||||
eo = toEigen(p1 - p0); //e2
|
||||
}
|
||||
// compute left and right cotangents
|
||||
double cotl = cotan(-el, eo); // -el -> angle between el and eo
|
||||
double cotr = cotan(er, eo);
|
||||
double cotl = cotan(-el, -eo);
|
||||
double cotr = cotan(-er, eo);
|
||||
// normalize edge vectors after cotangent computation
|
||||
el /= el.norm();
|
||||
er /= er.norm();
|
||||
|
||||
Loading…
x
Reference in New Issue
Block a user